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ABSTRACT 

We perform a set of fully self-consistent, dissipationless A^-body simulations to elucidate the dynamical re- 
sponse of thin galactic disks to bombardment by cold dark matter (CDM) substructure. Our method combines 
(1) cosmological simulations of the formation of Milky Way (MW)-sized CDM halos to derive the properties 
of substructure and (2) controlled numerical experiments of consecutive subhalo impacts onto an initially-thin, 
fully-formed MW type disk galaxy. The present study is the first to account for the evolution of satellite pop- 
ulations over cosmic time in such an investigation of disk structure. In contrast to what can be inferred from 
statistics of the z = Q surviving substructure, we find that accretions of massive subhalos onto the central regions 
of host halos, where the galactic disks reside, since z ^ I should be common. One host halo accretion history 
is used to initialize the controlled simulations of satellite-disk encounters. The specific merger history involves 
six dark matter substructures, with initial masses in the range 20% - 60% of the disk mass and of comparable 
size to the disk, crossing the central regions of their host host in the past ~ 8 Gyr. We show that these accre- 
tion events severely perturb the thin galactic disk and produce a wealth of distinctive dynamical signatures on 
its structure and kinematics. These include (1) considerable thickening and heating at all radii, with the disk 
thickness and velocity ellipsoid nearly doubling at the solar radius; (2) prominent flaring associated with an 
increase in disk thickness greater than a factor of 4 in the disk outskirts; (3) surface density excesses at large 
radii, beyond ^ 5 disk scale lengths, resembling those of observed antitruncated disks; (4) long-lived, lopsid- 
edness at levels similar to those measured in observational samples of disk galaxies; and (5) substantial tilting. 
The interaction with the most massive subhalo in the simulated accretion history drives the disk response while 
subsequent bombardment is much less efficient at disturbing the disk. We also explore a variety of disk and 
satellite properties that influence these responses. We conclude that substructure-disk encounters of the kind 
expected in the ACDM paradigm play a significant role in setting the structure of disk galaxies and driving 
galaxy evolution. 

Subject headings: cosmology: theory — dark matter — galaxies: formation galaxies: dynamics — galaxies: 
structure — methods: numerical 



1. INTRODUCTION 

Hierarchical models of cosmological structure formation, 
such as the c urrently favored cold dark matter (CDM ) 
paradigm (e.g.. lWhite & Rees|[T978t [Blumenflial et al.lll984 . 
generically predict substantial amounts of substructure in the 
form of small, dense, self-bound subhalos orbi ting within 
the virialized regions of larger host halos (e.g., Ghigna et aJJ 
■1998; Tormenet al. 1998; Moore et al. 1999; Klvpin et aT| 
1999bl: [Ghigna etal.ll2000l) . Observational probes' of sub- 
structure abundance thus constitute fundamental tests of 
the CDM model. Recently, a growing body of evi- 
dence has confirmed the hierarchical build-up of galaxy- 
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sized dark matter halos with the discovery of tidal streams 
and complex st ellar structures in the Milky Way (MW) 
(e.g..llbata et al .l 1994; Yannv et al. 2000; Ibata et al. 2001M 
Newberg et al. 2002; Maiewski et al. 2003; Martin et al 
2004; Martmez-Delgado et al. 2005; Grinmair & Dionatos 
2006; Belokurov et al. 2006), the Andro r neda galaxy (M31) 
jibata et al. 2001a; Fergus on etll.! 120021 l2005riKaUrai et al.l 
12006; Ibata et al. 2007), a nd beyon d the Local Group (e.g. , 
Malin & Hadlev 1997; Shang et al.l JT99l iPeng eTaLl 12003: 
Forbes et al. 2003; Pohlen et al. 200|). 

A significant fraction of observed galaxies have disk- 
dominant morphology with roughly 70% of Galaxy- 
sized dark matter halos hosting late -type systems (e.g., 
IWeinmann et al.l2006l;IChoi et al.l200"7h . Owing to the lack of 
a significant luminous component assoc iated with most sub- 
halos in Ga laxy-sized host halos (e.g., 'Klvpin et al."1999bt 
iMoore et"ar. 1999), information regarding the amount of sub- 
structure in these systems may be obtained via its gravitational 
influence on galaxy disks. Despite the small c ontribution of 
substructure to the total mass of the host (e.g., iGhigna et al.l 
2000), a considerable number of subhalos is expected within 
the virialized region of a Galaxy-sized CDM halo at any given 
epoch. If a large population of satellites exists, it may tidally 
disturb the host galactic disk, possibly leading to the imprint 
of distinctive dynamical signatures on its structure and kine- 
matics. 

Theoretical studies set within the CDM paradigm have 
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convincingly shown that the accretion of massive substruc- 
tures is com monplace during the formation of Galaxy-sized 
halos (e.g.. iLac ev & Co ld 119931: IZ entner & Bullockl 120031: 
iPurcell et al.H200 7: Stew art etaU boOS) and that typica l sub- 
halo orbits a re hig h ly eccentric ( e.g., Ghianaetal. 1998; 
iTormen et all 119981 : IZentner et all ©OSa; Benson 2005). 
These facts suggest that passages of massive satellites near 
the center of the host potential, where the galactic disks reside, 
should be common in CDM models. Such accretion events are 
expected to perturb the fragile circular orbits of disk stars by 
depositing large amounts of orbital energy into random stellar 
motions, gradually heating the disk and increasing its scale 
height. 

Yet, many disk galaxies are observed to be cold and 
thin, with average axial ratios of radial scale lengths 
to vertical sc ale heights, Rdlzd, in the range 4 - 5 
(e.g., de Griis 1998; Bizvaev & Mitronova 2002-. Kreg el etal.l 
1^02; Yoachim & Dalcanton 2006). In addition, recent stud- 
ies of edge-on disk galaxies using the Sloan Digital Sky 
Survey (SDSS) database have revealed a notable fraction of 
"super-thin" b ulgeless disks with much larger axial ratios 
dKautsch et al.ll2006i) . 

In the case of the MW, measurements of the disk 
scale height at the solar radius obtained using a vari- 
ety of methods, including star counts and mass mod- 
eling, indicate that the Galaxy comprises a thin, stel- 
lar disk with an exponential scale height of h. ~ 
300 ± 50 pc (e.g ., Kentetal. 1991; Dehnen & Binne^ 
1^98; Mendez & Guzman 1998; Larsen & Humphreys 2003t 
Widro w & Dub inski 2005; Juric et al. 2008). Furthermore, 



the age-velocity dispersion relation of disk stars in the 
solar neighborhood suggests that a significant fraction of 
the thin disk of the MW was in place by z ^ 1 (e.g. , 
Wvse"2001^, 'Guillen & G mietillOOOt iNordstrom et al.ll2004t 
Seabroke & Gilmore 200f~ The existence of such an old, 
thin stellar disk, as established by the age distribution of disk 
stars, may imply an absence of satellite accretion events over 
the past ^ 8 Gyr. Such an extended period of quiescent dy- 
namical evolution is difficult to reconcile with the hierarchi- 
cal assembly of structure prescribed by the ACDM cosmo- 
logical model. Although hydrodynamical simulations of disk 
galaxy formation have offered some insights into accommo- 
datin g observational fac ts that challenge the CDM paradigm 
(e.g., Abadi et ani2003h . the detailed dynamical response of 
galactic disks to halo substructure in a cosmological context 
remains poorly understood. 

Significant theore tical effort, in c luding both semi- 
analy ti c modeling dToth & Ostrikeii 119921 iBenson et alj 
120041; iHop kins et al.l 20081) and nurneric al simulations 
dOuinn & G oodman ligset IOuinn et"an [199 3; Walker et al] 
1996"; 'Huang & CarTber? '1997; 'Sellwood et al. 199a 
Velazquez & White 1999;' Font et al. 2001; Ardi et al. 200l 
Gauthier et al. 2006; Havashi & Chiba 2006; Rea d et all 
200& .Villalobos & Helm! 2008; Purcell et al. 2008) has been 



devoted to quantifying the resilience of galactic disks to 
infalling satellites. Although valuable in several respects, 
these earlier investigations could not capture fully the amount 
of global dynamical evolution induced in thin disk galaxies 
by substructure in the context of the ACDM model. While a 
detailed comparison to previous work is presented in § |6l we 
mention at the outset that each of the aforementioned numer- 
ical studies suffered from at least one critical shortcoming by 
either: (1) not being fully self-consistent, modeling various 
components of the primary disk galaxy and/or the satellites 



as rigid potentials, a choice which leads to overestimating 
the damage caused to the disk; (2) focusing on experiments 
with infalling systems on nearly circular orbits that are poor 
approximations of the highly eccentric orbits typical of CDM 
substructure; (3) adopting galactic disks that are much thicker 
compared to typical thin disks including the old, thin stellar 
disk of the MW; (4) modeling the compact, baryonic cores of 
accreting systems exclusively and neglecting the more diffuse 
and extended dark matter component; and (5) considering the 
encounters of individual satellites with galactic disks, despite 
the CDM expectations of numerous accretion events over the 
history of a galaxy. 

Reg arding the final po int, notable exceptio ns were the stud- 
ies bv lFont et all d2001h and lGauthier et alj ([2006) which ex- 
amined the dynamical evolution of a stellar disk in the pres- 
ence of a large ensemble of dark matter subhalos. Both con- 
tributions reported negligible tidal effects on the global struc- 
ture of the disks. However, these investigations had the draw- 
back of adopting the z = surviving substructure present in 
a Galaxy-sized CDM halo, and thus not accounting for past 
encounters of systems with the galactic disk during the evo- 
lution of the satellite populations. This point is critical be- 
cause as subhalos on highly eccentric orbits continuously 
lose mass, the number of massive satellites with small or- 
bital pericenters that are most capable of severely perturbing 
the disk decli nes with redshift so that few would be present 
at z = (e.g., 'Zentner & Bullockl l200l lKravtsov""etaLll2004l; 
iGao et al. 2004). Establishing the role of halo substructure in 
shaping the fine structure of disk galaxies clearly requires a 
more realistic treatment of its evolution over cosmic time. 

Our aim is to improve upon these shortcomings and extend 
consideration to the rich structure of perturbed galactic disks. 
This paper is the second in a series elucidating the effects 
of halo substructure on thin disk galaxies in the context of 
the prevailing CDM paradigm. Kazantzidis et al. (2008, here- 
after Paper I) focused on the generic morphological signatures 
induced in the disk by a typical ACDM-motivated satellite 
accretion history, while the present paper discusses the dy- 
namical response of the galactic disk subject to bombardment 
by the same population of dark matter subhalos. We imple- 
ment a two-step strategy in an effort to overcome some of the 
drawbacks of past studies. First, we analyze dissipationless, 
cosmological simulations of the formation of Galaxy-sized 
CDM halos to derive the accretion histories and properties 
of substructure populations. This information is subsequently 
used to seed collisionless, controlled numerical experiments 
of consecutive satellite impacts onto A^-body realizations of 
fully-formed, thin disk galaxies. Given the outstanding issues 
regarding di sk galaxy formation in CDM cosmogonies (e.g., 
iMaver et al. 2008) as well as the inadequacies of the current 
generation of cosmological simulations to resolve all dynam- 
ical scales and physical processes relevant to satellite-disk in- 
teractions, this strategy is most appropriate. As in Paper I, we 
model the infalling systems as pure dark matter subhalos and 
focus exclusively on the evolution of the stellar material in the 
disk itself. 

Our contribution improves upon earlier studies in several 
important respects. First and foremost, we examine the re- 
sponse of galactic disks to subhalo populations that are truly 
representative of those accreted and possibly destroyed in the 
past. As such, we mitigate the biases in the incidence of ac- 
cretion events and properties of infalling satellites induced 
by considering only the z = substructure of a CDM halo. 
Specifically, we extract merger histories of host halos since 
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z 1 and study the ramifications of such accretion events for 
disk structure. As we illustrate below, our methodology re- 
sults in a substantially larger number of potentially damaging 
subhalo-disk encounters than previously considered. Second, 
we construct self-consistent satellite models whose properties 
are culled directly from the same cosmological simulations 
of Galaxy-sized CDM halos. This obviates the need for ar- 
bitrary assumptions regarding the numbers, masses, internal 
structures, orbital parameters, and accretion times of the in- 
falling subhalos. 

Lastly, we employ multi-component disk galaxy models 
that are derived from explicit distribution functions, are mo- 
tivated by the ACDM paradigm, and are flexible enough to 
permit detailed modeling of actual galaxies such as the MW 
and M3 1 by incorporating a wide range of observational con- 
straints. These properties in synergy with the adopted high 
mass and force resolution enable us to construct realistic, 
equilibrium A^-body models of thin disk galaxies and study 
their dynamical response to encounters with satellites. For 
our primary simulation set, we use a best-fit model for the 
present-day structure of the MW and employ an initial disk 
scale height of Zd = 400 pc which is consistent with that of 
the old, thin stellar disk of the Galaxy. Although we utilize 
a model for the MW, we do emphasize that our simulation 
campaign is neither designed to follow the evolution of nor to 
draw specific conclusions about the Galaxy or any other par- 
ticular system. Given the complex interplay of effects (e.g., 
gas cooling, star formation, chemical evolution) relevant to 
the formation and evolution of spiral galaxies, our collision- 
less simulations aim for generic features of the evolution of a 
thin galactic disk subject to bombardment by CDM substruc- 
ture 

Our work establishes that the types of merger histories ex- 
pected in ACDM can substantially perturb the structure of a 
cold, stellar disk. We demonstrate that cosmological halo as- 
sembly via multiple accretion events is capable of generating 
a wealth of distinctive dynamical signatures in the structural 
and kinematic properties of disk stars. These include pro- 
nounced thickening and heating, prominent flaring and tilting, 
surface density excesses that develop in the outskirts similar 
to those of observed antitruncated disks, and lopsidedness at 
levels similar to those measured in observed galaxies. Our 
findings suggest that details of the galaxy assembly process 
may be imprinted on the dynamics of stellar populations and 
corroborate the concept of accretion-induced galaxy evolu- 
tion. We also show that the global dynamical response of a 
galactic disk to interactions with satellites depends sensitively 
on a variety of parameters including the initial disk thickness, 
the presence of a bulge component in the primary disk, the 
internal density distribution of the infalling systems, and the 
relative orientation of disk and satellite angular momenta. 

The remainder of this paper is organized as follows. In §|2] 
we describe the methods employed in this study. § |3] con- 
tains the results regarding the dynamical signatures induced 
in a thin galactic disk by a typical ACDM accretion history. 
In § m we study various factors that could influence the re- 
sponse of disk galaxies to satellite accretion events. Implica- 
tions and extensions of our findings along with a comparison 
to previous work are presented in §|5]and §|6] Finally, in §|2l 
we summarize our main results and conclusions. Through- 
out this work we use the terms "satellites," "subhalos," and 
"substructures" interchangeably to indicate the distinct, grav- 
itationally bound objects that we use as the basis for the con- 
trolled satellite-disk encounter simulations. 



2. METHODS 

A thorough description of the adopted methodology is 
given in Paper I. For completeness, we provide a summary of 
our approach here and refer the reader to Paper I for details. 

2.1. Hierarchical Cosmological Simulations 

We analyze high-resolution, collisionless cosmological 
simulations of the formation of four Galaxy-sized ha- 
los in a flat ACDM cosmological model with parameters 
(ri„,fiA,nfc,/j,CT8) = (0.3,0.7,0.043,0.7,0.9). The simula- 
tions were performe d with the Adaptive R efinement Tree 
(ART) N-hoAy (Ki'a vtsov et all [19971 iKravTs ov 1999). The 
host halos that we considered come from two different sim- 
ulations. Halos "Gi," "G2," and "G3" all formed within a 
cubic box of length 25 /!"'Mpc on a side, while halo "G4" 
formed in a cubic volume of 20 /i~'Mpc on a side. The 
mass and force resolution of these simulations as well as var- 
ious properties of halos G1-G4 and their substructure were 



discussed in Paper I and in previous literature ( Klvni n et al 



2001 ; [Kravtsov et al. 2004; Zentneret al. 2005b; Pradaetal 
200d; iGnedin & Ki-avtsov .2006) . We identify halos (out- 



side hosts G1-G4) and substructur es using a va riant of the 
Bound Density Maxima algorithm (Klypin et alJ[l999 a) and 
we have constructed detailed accretion histories for each host 
halo and orbital tracks for all subhalos. The specifics of this 
analysis can be foun d in Paper I and are based largely on 
iKravtsov et"an ( 120041) . 

Our conventions are as follows. We adopt a mean over- 
density of 337 to define the virial radius, rvh, and the cor- 
responding virial mass, Mvir, for each host halo at z = 0. 
Halos Gi through G4 have r™ ~ (234,215,216,230) /j-'kpc 
andMvir -(1.5,1.1, 1.1, 1.4) x lO'^ /z-'Mq, respectively. All 
four of these halos accrete only a small fraction of their fi- 
nal mass and experience no major mergers at z < 1, and 
are thus likely to host a disk galaxy. Moreover, their accre- 
tion histories are typ ical of systems in this mass range (e.g., 
IWechsler et al.l l20^h . We have chosen Gi as our fiducial 
case for the satellite-disk interaction experiments described 
in §1231 

2.2. Interactions Between Substructures and Disks in CDM 

We investigate the dynamical response of thin galactic disks 
to interactions with CDM substructure incorporating for the 
first time a model that accounts for its evo lution over cos- 
mic time. Previous related studies (e.g.. Font et akl 120011 : 
lGauthier"er al. 2006) only utilized subhalo populations at z = 
0, rather than a complete merger history, to seed simulations 
of satellite-disk encounters in a cosmological context. This 
procedure has the drawback of eliminating from considera- 
tion those massive satellites that, prior to z = 0, cross through 
the central regions of their hosts, where the galactic disks re- 
side. Such systems can potentially produce strong tidal effects 
on the disk, but are unlikely to constitute effective perturbers 
at z = as they suffer substantial mass loss (or even become 
disrupted) during their orbital evolution p recisely b ecause of 
their forays into the central halo (e.g . . IZentner & Bulloc 
2003; Ki-avtsov et al.' I2004t iGaoetalJ 120041: IZentner et ' 
2005a; Benson 2005). Accounting for the impact of these rel- 
atively short-lived objects on the global dynamical response 
of galactic disks is the major improvement we introduce in 
the present study. 

Figure [T] serves as a dramatic illustration of this element 
of our modeUng. This figure is a scatter plot of mass ver- 
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Fig. 1. — A scatter plot of mass versus pericentric distance for sublialos 
identified in four Galaxy-sized dark matter halos formed in the ACDM cos- 
mology. Results are shown after scahng the virial quantities of all host halos 
to the corresponding values of the primary disk gala xy model used in the con- 
trolled satellite-disk encounter simulations of § 12.41 Subhalo masses and peri- 
centers are assigned according to the description in the text and are presented 
in units of the mass, Mjisk = 3.53 X IO'^Mq, and scale length, = 2.82 kpc, 
respectively, of the disk in the same galaxy model. Filled symbols refer to 
systems that cross within a (scaled) infall radius of ri„f = 50 kpc from their 
host halo center since a redshift z = 1 . Symbols with crosses correspond to the 
specific subhalos in host halo G; used to seed the controlled simulations of 
satellite-disk interactions. Open symbols refer to the z = surviving substruc- 
tures. Dotted lines mark the so-called "danger zone" with Msub ^ 0.2Mdjsk 
and Tpeii < 20 kpc corresponding to infalling subhalos that are capable of 
substantially perturbing the disk. Close encounters between massive sub- 
structures and galactic disks since z = 1 should be common occurrences in 
ACDM. In conttast, very few satellites in present-day subhalo populations 
are likely to have a significant dynamical impact on the disk structure. 

sus pericentric distance for two different substructure popu- 
lations within host halos G1-G4. The masses and pericen- 
ters of all subhalos have been scaled to the mass, Mdisk = 
3.53 X IO'^Mq, and radial scale length, Rd = 2.82 kpc, of 
the disk in the primary galaxy model used in the controlled 
satellite-disk encounter simulations of § 12.41 For the purposes 
of this presentation, we have also scaled the virial quantities 
(Mvir and ryir) of all four Galaxy-sized host halos to the total 
mass. Ml, = 7.35 x 1O"M0, and tidal radius, Ri, = 244.5 kpc, 
of the parent dark matter halo in the same galaxy model. The 
full list of parameters pertaining to the primary disk galaxy 
will be listed in § 12331 

The first substructure population of Figure[T]comprises sys- 
tems that cross within a (scaled) infall radius of = 50 kpc 
from their host halo center since a redshift z = 1 . This se- 
lection is fixed empirically to identify orbiting satellites that 
approach the central regions of the host potential and are thus 
likely to have a significant dynamical impact on the disk struc- 
ture (Paper I). We assign masses to the satellites of this group 
at the simulation output time nearest to the first inward cross- 
ing of rinf . As we discuss below, we define this to be the epoch 
that our controlled simulations initiate. The corresponding 
pericen ters are computed f rom the orbit of a test particle in 
a static iNavarro et al.l (119961 hereafter NEW) potential whose 
properties match those of the host CDM halo at the time of 



Hnf. We note that a single distinct object of this population 
may be recorded multiple times as one subhalo may undergo 
several passes through the central regions of its host with dif- 
ferent masses and pericenters. Many of these satellites suffer 
substantial mass loss or even become tidally disrupted prior 
to z = 0. 

The second subhalo population consists of the z = Q surviv- 
ing substructures. Their pericenters are also estimates based 
on the orbit of a test particle in a static NEW potential whose 
properties match those of the host CDM halo at z = 0. The dot- 
ted line in Eigure [T] encloses an area in the Msub-'"peri plane 
corresponding to satellites more massive than 0.2Mdisk with 
pericenters of rped < 20 kpc (rperi < 7/?^/). We refer to this 
area as the "danger zone". Subhalos within this area should 
be effective perturbers, but we intend this as a rough criterion 
to aid in illustrating our point. 

Eigure [T] demonstrates that the z = substructure popula- 
tions contain very few massive systems on potentially dam- 
aging orbits. In fact, in all four of our host halos only one 
satellite can be identified inside the danger zone at z = 0. 
As a consequence, the present-day substructure in a CDM 
halo likely plays only a minor role in driving the dynamical 
evolution of galactic disks. This co n clusi on is supported by 
iFont et all (12001) and Ga uthier et alJ fcOOd) who did consider 
the z = subhalo populations and reported negligible effects 
on the global structure of the disk. 

On the other hand, the danger zone contains numerous 
satellites that have penetrated deeply into their host halo since 
z = 1: on average, ^ 5 systems more massive than 0.2Mdisk 
pass close to the center of their hosts with rperi S 20 kpc in the 
past ^ 8 Gyr We also stress that three of the host halos have 
accreted at least one satellite more massive than the disk itself 
since z = 1 . This finding is corroborated by analysis of simula- 
tions with much better statistics (iStewart et alJ 2008). Overall, 
the results in Eigure[T]indicate that close encounters between 
massive subhalos and galactic disks since z = 1 should be com- 
mon occurrences in the ACDM cosmological model. It is thus 
important to quantify and account for such interactions when 
the goal is to investigate the cumulative dynamical effects of 
halo substructure on galactic disks. 

2.3. Initial Conditions for Satellite-Disk Encounters 

Initializing the controlled experiments of subhalo-disk en- 
counters involves: (1) identifying relevant substructures in the 
cosmological simulations and recording their properties (mass 
functions, internal structures, orbital parameters, and accre- 
tion times); and (2) constructing A^-body realizations of disk 
galaxy and satellite models. In what follows, we describe each 
of these steps. 

2.3.1. Subhalo Selection Criteria 

In practice, it is not feasible to simulate the encounters of 
all orbiting satellites present in host halo Gi with the galactic 
disk at sufficient numerical resolution. This, however, is not 
needed as the dynamical effects of substructure on galactic 
disks are expected to be dominated by the few most massive 
subhalos that penetrate deeply into the central regions of their 
hosts. To limit computational expense as well as to ensure 
that only relevant encounters will be considered, we impose 
two criteria for selecting target satellites for re-simulation (see 
also Paper I). 

We limit our search to satellites that approach the central 
regions of halo Gi with (scaled) pericenters of rperi ^ 20 kpc 
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since z = 1 . This choice is motivated by the fact that sub- 
halos with small orbital pericenters are expected to substan- 
tially perturb the galactic disk. We initiate our controlled 
re-simulations at the epoch when each selected subhalo first 
crossed a (scaled) infall radius of ri^f = 50 kpc from the host 
halo center. In what follows and unless otherwise explicitly 
stated all satellite properties (obviously except for those at 
z = 0) are recorded at the simulation output closest to the first 
inward crossing of ri„f. Likewise, we select only subhalos 
that are a significant fraction of the disk mass, with (scaled) 
masses of M^ub ^ 0.2Mdisk, as these will have the largest im- 
pact on the disk structure. 

Figure [T] shows that, on average, ^ 5 systems meet these 
two criteria for each host halo. Although at least one sub- 
halo more massive than the disk in our controlled simulations 
is expected to cr oss the central reg ions of their hosts since 
z ~ 1 (Figure [1] iStewart et al.l 12008), we explicitly ignore 
these interactions with Mjub ^ Mdisk- Our goal is to investi- 
gate the tidal effects of substructure on a well-preserved disk 
galax y while such encou nters may destroy the disk upon im- 
pact ( Purcell et al.]|2008l) . 

The aforementioned criteria result in a merger history of 
six accretion events, which we denote S1-S6, for subsequent 
re-simulation over a « 8 Gyr period. These substructures cor- 
respond to the filled symbols with crosses seen in Figure[T]and 
their properties are summarized in Table [T] below. A critical 
reader may note that the internal properties and orbital param- 
eters of subhalos in collisionless cosmological simulations, as 
those of the present study, may differ substantially from those 
in simulations with gasdynamics. Indeed, in hydrodynami- 
cal simulations of disk galaxy formation, the average satel- 
lite crossing the central regions of its host experiences addi- 
tional mass loss due to disk shocking and both circularization 
of its orbit and decrease of its orbital inclination via dynam- 
ical friction aga inst the di sk (e.g. , Ouinn & Goodman 1986; 
iPenarrubia et al. 2002; Meza et al.il2005h . These considera- 
tions do not apply to subhalos S1-S6 as these systems were 
identified in the cosmological simulation and their properties 
were recorded before crossing the center of their host halo for 
the first time. 

Table [1] contains the main orbital and structural properties 
of cosmological satellites S1-S6. The parameters of each sub- 
halo in this table are scaled so that they correspond to the 
same fractions of virial quantities (Mrvir, r^m, and Vyir) in both 
host halo Gi and in the primary galaxy model. Columns 
(2)-(9) list properties measured in the cosmological simula- 
tion. Columns (10)-(12) present structural parameters com- 
puted directly from the A^-body realizations of satellite models 
constructed for the controlled subhalo-disk encounter experi- 
ments. 

The properties of satellites S1-S6 have been discussed ex- 
tensively in Paper I. Nonetheless, a few parameters are worthy 
of explicit note. First, the masses of these six substructures 
correspond to the upper limit of the mass function of observed 
satellites in the Local Group. For reference, the mas s of the 
Large Magellanic Cloud (e.g.. lSchommer et al.lll992l) is sim- 
ilar to that of our most massive subhalo S2. Moreover, satel- 
lites S 1 -S6 are spatially extended with their tidal radii, rtid, en- 
compassing a significant fraction of the disk itself (rtid ^ TRd)- 
Because of that the entire disk will be subject to potential fluc- 
tuations, although of different magnitude in different places, 
during the encounters with these subhalos. Thus, the energy 
imparted by typical cosmological substructures will likely not 
be deposited locally at the point of impact, as it was assumed 




Fig. 2. — Upper panel: Spherically-averaged density profiles, p(r), for 
representative cosmological satellites S 1 and S2 as a function of radius in 
units of the tidal radius of each system, rtij. Stars correspond to the sub- 
halo profiles extracted directly from the cosmological simulation of host halo 
Gi. Solid lines present fits to the density structure using a multi-parameter 
(«,/3,7) density law, while dotted lines show coiresponding NFW fits. The 
solid lines include an exponential cutoff at r > r^d. All curves are plotted 
from the adopted force resolution (2es„b = 300 pc) outward and densities are 
normalized to the mean density within the tidal radius of each satellite. For 
clarity, the density profiles corresponding to the lower curves are vertically 
shifted downward by a factor of 0.05. The (a, l3, 7) density law provides an 
accurate description of the internal structure of cosmological subhalos at all 
radii, while the NFW functional form substantially underestimates subhalo 
densities in the innermost regions. Bottom panel: Residuals for the density 
profile fits, i^p I p = Psub ~Pfl(/Psub> where ftuj, is the true subhalo density 
computed in the cosmological simulation and pu, denotes the fitted density. 



bv iToth & Ostrike^ (119921) , but rather globally across the en- 
tire disk. Lastly, we note that most subhalos S1-S6 are on 
highly eccentric orbits and that the simulated accretion his- 
tory includes nearly polar (S1,S2), prograde (S3,S5), and ret- 
rograde (S4,S6) encounters. 

2.3.2. Constructing N-body Realizations of Satellites 

Our controlled experiments employ satellite models which 
are constructed to match the internal structure of cosmolog- 
ical subhalos S1-S6. The density profiles of these systems 
are extracted at the simulation output closest to each sub- 
halo's first inward crossing of ^nf. We model all cosmolog- 
ical satellites with the general density profile law (IZhaoll996t 
.Kravtsov et al. 1998) 



p{r) = 



{r < nid) , (1) 



where sets the normalization of the density profile and r^d 
is the limiting radius of the satellite imposed by the host halo 
tidal field. The exponents 7 and (3 denote the asymptotic in- 
ner and outer slopes of the profile, respectively. Exponent a 
parametrizes the transition between the inner and outer pro- 
files, with larger values of a corresponding to sharper tran- 
sitions. Lastly, at the scale radius rj, the logarithmic slope 
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TABLE 1 

Parameters of the Satellite Models 



Model 


Z-dCC 




rtid/Rd 




''apo/''peri 


(°) 


(kms"') 




(a, A 7) 


(kpc) 


c 


(1) 


(2) 


(3) 


(4) 


(5) 


(6) 


(V) 


(8) 


(9) 


(10) 


(11) 


(12) 


SI 


0.96 


0.33 


8.8 


2.6 


6.8 


93.3 


72.9 


0.41 


(0.12,3.85.1) 


2.1 


11.6 


S2 


0.89 


0.57 


7.6 


2.6 


6.0 


86.6 


70.9 


0.46 


(0.21,4.02,1) 


2.6 


8.1 


S3 


0.54 


0.42 


8.2 


6.2 


3.2 


45.1 


158.2 


0.72 


(0.38,3.72,1) 


2.2 


10.6 


S4 


0.32 


0.45 


7.0 


0.5 


20.6 


117.7 


19.6 


0.16 


(0.25,4.18,1) 


1.6 


12.5 


S5 


0.20 


0.22 


9.7 


3.7 


9.3 


59.9 


171.3 


0.35 


(0.16,3.94,1) 


1.8 


15.2 


S6 


0.11 


0.21 


8.2 


1.1 


19.6 


144.5 


89.6 


0.17 


(0.29,4.09,1) 


1.3 


18.3 



is the average of the inner and outer slopes, dlnp(r)/dlnr = 

-(7+/3)/2. 

In practice, all subhalos are described well by Eq. ([T]i with 
7 = 1, so we fix this value and fit their structure by varying 
the other parameters. Moreover, we truncate subhalo profiles 
beyond the tidal radius with an exponential law 



p(r) = p(rtid) 



Aid 



exp 



^ decay 



(r > Aid) , (2) 



Notes. — Columns (2)-(9) record satellite properties at the epoch closest to when each subhalo first crossed within a (scaled) infall radius of rinf = 50 kpc 
from the center of host halo Gi . Colum ns (1 0)-(12) list structural parameters computed from the A'-body realizations of satellite models used in the controlled 
subhalo-disk encounter simulations of § 12.41 All entries are listed after scaling the virial quantities of halo Gj to the corresponding values of the parent halo in 
the fiducial disk galaxy model used in the controlled experiments. Note that values of orbital circularities and apocenter-to-pericenter ratios may deviate from 
those measured directly in the controlled simulations because the former are estimated from the satellite orbits in the potential of host halo Gi . Col. (1): Satellite 
model. Col. (2): Redshift at which the properties of cosmological satellites were recorded. Col. (3): Bound satellite mass in units of the mass of the disk, 
Mdi,sk = 3.53 X 1O"'M0, in the fiducial galaxy model used in the controlled encounter simulations. Col. (4): Satellite tidal radius in units of the radial scale 
length of the fiducial disk, Rj = 2.82 kpc. Col. (5): Pericenter of the satellite orbit in units ofRj. Col. (6): Satellite orbital apocenter-to-pericenter ratio. Col. (7): 
Angle between the initial angular momenta of the satellite and the disk in degrees. This angle is defined so that 0° < 8 < 90° corresponds to a prograde orbit 
and 90° < 6 < 1 80° corresponds to a retrograde orbit. The cases where 6 = 0°, 6 = 90°, and 8 = 180° represent a coplanar prograde, a polar, and a coplanar 
retrograde orbit, respectively. Col. (8): Satellite 3D orbital velocity in kms"'. Col. (9): Circularity of the orbit. This parameter is defined as tj = //7(;irc(£). 
where / is the satellite angular momentum and Jchc(E) is the corresponding angular momentum for a circular orbit of the same energy E. Col. (10): Intermediate, 
outer, and inner slopes of the satellite density profile. Inner and outer slopes correspond to the asymptotic values. Col. (11): Scale radius of the satellite density 
profile in kpc. This radius is defined as the distance where the logarithmic slope is the average of the inner and outer slopes, dlnp(r)/dlnr = -(7 + f3)/2. Col. 
(12): Satellite concentration defined as c ~ rty/cs. 

ical subhalos. This is important because, as we illustrate in 
§ 14.31 the dynamical response of galactic disks to encounters 
with satellites depends sensitively upon the precise density 
distribution of the infalling systems. 

We list fit parameters for each subhalo in Table[TJ however, 
we note two general features. First, in all cases the asymptotic 
outer slopes (3 are found to be much steeper than 3. Best- 
fit values vary from /3 = 3.7 to /? = 4.2 suggesting that these 
satellites are better approxi mated by the Hernquist ( 1990) pro- 
file at large radii (see also iGhigna et al.i,1998.) . Second, the 
intermediate slopes a range from a = 0.1 -0.4 indicating a 
much more gradual transition between inner and outer asymp- 
totic power-laws compared to that of the NEW profile with 
a= I. The steep outer profiles and variety in structural pa- 
rameters is, at least in part, due to the strong dynamical evolu- 
tion that subhalos e xperience within the i r host potential (e.g. , 
iMoore et al .l 119961; iK lvpin et al."1999a': 'Havashi et al."2003l: 
Kazantzidi s et alj l200 4b; Kravtsov et al. 2004; Maver et al] 
2007). 

A^-body realizations of satellites are constructed from a dis- 
tribution function (DF) that self-consistently reproduces the 
density structures of selected subhalos S1-S6. Because sub- 
structures in cosmological simulations a re nearly spherical 
in both configuration and velocity space (M oore et alj|2004l ; 
Kazantzidis et al. 2006; Kuhlen et al. 2007), we assume that 
the DE depends only upon the energy per un it mass and cal - 
culate it through an Abel transform ( Kazantzidis et alj2004al) . 

Lastly, each satellite is represented by A^sub =10^ particles. 
The gravitational softening length is set to e^ub = 150 pc which 
allows us to resolve the structure of subhalos to ~ 1 % of their 
tidal radii. The adopted mass and force resolution are ade- 
quate to resol ve mass loss processes ex perienced by orbiting 
satellite halos ( [Kazantzidis et al.ll2004bl) . 



where k is fixed by the requirement that dlnp(r)/dlnr is con- 
tinuous at rtid- This procedure is necessary because sharp trun- 
cations result in subhalo models that are not in equilibrium 
dKazantzidis et al J 120 04a'), but it results in additional bound 
mass beyond rtid. The precise amount of additional mass de- 
pends upon the model parameters, but is roughly ^ 2% of the 
bound mass for each subhalo that we consider 

Figure |2] demonstrates the efficacy of our procedure for 
modeling the internal structure of cosmological subhalos. 
This figure presents spherically-averaged density profiles, 
p(r), of two representative cosmological satellites along with 
two different fits to their density distributions and the associ- 
ated residuals, Ap/ p. The first is a fit to the multi -parameter 
(a,/?, 7) functional form introduced above [Eq. ([T]i], while 
the second is a fit to the NEW density profile with (a,/?, 7) = 
(1,3,1). 

Figure |2] shows that (a,/?, 7) models provide accurate rep- 
resentations of the density structures of cosmological subha- 
los, while the NEW profile much less so, as one might expect 
because it has less parameter freedom. Indeed, the density 
residuals with respect to Eq. ([T]) are < 10% for both satellites 
over the entire range of radii, while the residuals to NEW fits 
are much larger reaching ^ 50% in the innermost parts of the 
profile. We note that these findings were confirmed in the re- 
maining subhalos S3-S6. Overall, Figure |2] demonstrates that 
the NEW law underestimates the normalization of the density 
profile and consequently the concentration of these cosmolog- 
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2.3.3. Disk Galaxy Models 

We construct A^-body realizations of multi-component pri- 
mary d isk galaxies using the method of IWidrow & Dubinskil 

(120051) as describ ed in detail in Paper I. The 

IWidrow & DubinskH (I2005h models are derived from 
explicit DFs and thus represent axisymmetric, equilibrium 
solutions to the coupled coUisionless Boltzmann and Poisson 
equations. Because of that, they are ideal for studying 
complex dynamical processes associated with the intrinsic 
fragility of galactic disks such as gravitational interactions 
with infalling satellites. The galaxy models consist of a n 
exponential stellar disk, a Hernquist bulge ( Hernquist 1990h . 
and an NFW dark matter halo, and are characterized by 15 
free parameters that may be tuned to fit a wide range of 
observational data for actual galaxies. 

For the majority of satellite-disk encounter simulations, 
we use the specific parame ter choices of model "MWb" in 
IWidrow & Dubinskil (l2005 h. which satisfies a broad range of 
observational constraints on the MW galaxy. The stellar disk 
has a mass of Mdisk = 3.53 x IO'^Mq, a radial scale length 
of Rd = 2.82 kpc, and a sech^ scale height of Zd = 400 pc. 
The latter is consistent wit h that inferred for the old, thin stel- 
lar disk of the MW (e.g., Kent et al. 1991; Dehnen & Binnev' 
1998; Men dez & Guzman 1998; Larsen & Humphreys 2003; 
Juric et aP l2008l) . It is reasonable to initialize the galac- 
tic disk with suc h thickn ess as observational e vidence (e.g., 
lOuillen & Garne tt 2000; Nordstrom et al.l2004l) indicates that 
the scale height of the thin disk of the MW has not changed 
significantly over the period modeled in this study (z < 1). 
The bulge has a mass and a scale radius of Mi, = 1.18 x 
10"*Mq and a/, = 0.88 kpc, respectively. The dark matter 
halo has a tidal radius of Rh = 244.5 kpc, a mass of Mh = 
7.35 X IO'^Mq, and a scale radius of n, = 8.82 kpc. The to- 
tal circular velocity of the galaxy model at the solar radius, 
Rq ~ 8 kpc, is VciRQ) = 234.1 kms"' and the Toomre disk 
stability parameter is equal to Q = 2.2 at R = 2.5Rd- We note 
that direct numerical simulations of the evolution of model 
MWb in isolation confirm its stability against bar formation 
for 10 Gyr. Therefore, any significant bar growth and asso- 
ciated disk evolution identified during the satellite-disk en- 
counter experiments should be the result of subhalo bombard- 
ment. In what follows, we refer to this galaxy model as "Dl." 

In addition to our primary simulation set, we also address 
the dependence of disk response to encounters with infalling 
subhalos upon initial disk thickness and the presence of a 
bulge. Thicker disks have larger vertical velocity dispersions 
so we might expect them to be more robust to satellite ac- 
cretion events. In addition, a central bulge component may 
act to reduce the amount of damage done to the structure of 
the inner disk by the infalling subhalos. Correspondingly, we 
initialize two additional disk galaxy models. 

The first modified galaxy model was constructed with the 
same parameter set as Dl, except that it was realized with a 
scale height that was larger by a factor of 2.5 {zd = 1 kpc). 
We refer to this "thick" disk galaxy model as "D2." Except 
from disk thickness, all of the other gross properties of the 
three galactic components of model D2 are within a few per- 
cent of the corresponding values for Dl. The second modi- 
fied galaxy model is the same as Dl, but constructed without 
a bulge component. We refer to this "bulgeless" disk galaxy 
model as "D3." We stress that a bulgeless ver sion of Dl con- 
structed using the 'W idrow & Dubinskil ( 120051) method differs 
significantly from Dl. This is because the DF of the compos- 



ite galaxy model is related to the individual density distribu- 
tions and DFs of each galactic component in a non-trivial way. 
To mitigate such differences which make model comparison 
cumbersome, we have chosen to realize disk model D3 by 
adiabatically evaporating the bulge component from galaxy 
model Dl. The timescale of evaporation was set to 500 Myr 
to minimize the response of the disk and halo to the removal 
of the bulge. During the evaporation process, while the main 
properties of the inner disk and halo (e.g., density profiles, ve- 
locity ellipsoids) have changed in response to the decrease of 
the central potential, the disk thickness has remained largely 
unmodified. It is also worth mentioning that when evolved 
in isolation model D3 develops a bar inside ~ 4 kpc at time 
t - 2 Gyr 

For each disk galaxy model, we generated an A^-body re- 
alization containing Nd = 10^ particles in the disk and Ni, = 
2x10^ particles in the dark matter halo. In experiments with 
models Dl and D2, bulges were represented with A';, = 5 x 10^ 
particles. The gravitational softening lengths for the three 
components were set to = 50 pc, e/, = 100 pc, and e/, = 50 pc, 
respectively. Mass and force resolution are sufficient to re- 
solve the vertical structure of our galactic disks as well as 
minimize artificial heating of the disk particles through in- 
teractions with the much more massive halo particles. Lastly, 
in all cases we oriented the primary galaxy models such that 
the disk and host halo Gi angular momenta were aligned (e.g., 
iLibeskind etani2007l) . 

2.4. Satellite-Disk Encounter Simulations 

All controlled simulations of satellite-disk interactions 
were carried out with the PKDGRA V multi-stepping, paral- 
lel, tree A^-body code dStadell 120011) as described in Paper I. 
We treated impacts of cosmological subhalos S1-S6 onto the 
disk as a sequence of encounters. Recall that the internal prop- 
erties and orbital parameters of these satellites were scaled so 
that they correspond to the same fractions of virial quantities 
in both host halo Gi and in the primary disk galaxy model. 
Starting with subhalo SI we included subsequent systems at 
the epoch when they were recorded in the cosmological sim- 
ulation (Table [U. This procedure is described in more de- 
tail in Paper I. We note here that the sequence of accretion 
events is such that S 1 and S2 undergo simultaneous interac- 
tions with the disk. We have conducted an additional exper- 
iment in which satellite S2 was introduced only after the in- 
teraction with S 1 was completed and confirmed that the disk 
dynamical response is similar in the two cases. 

Due to their considerable mass and size, the infalling sub- 
halos were not introduced in the simulations directly; instead, 
they were grown adiabatically in their orbits. This method 
ensures that the disk does not suffer substantial perturbations 
due to the sudden presence of the satellite and change in po- 
tential at its vicinity. Specifically, each encounter simulation 
was performed using the following procedure: (i) Insert a 
massless particle realization of each satellite at the distance 
at which it was recorded in the cosmological simulation of 
host halo Gi. (ii) Increase the mass of this distribution to 
its final value linearly over a timescale that ranges between 
~ 150 and 400 Myr depending on subhalo mass. During 
this period, the satellite remained rigid and its particles are 
fixed in place, while the dark and baryonic components of the 
primary disk galaxy were allowed to achieve equilibrium with 
the subhalo as its mass grows, (iii) Initialize the "live" A^-body 
satellite model by setting its internal kinematics according to 
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the method described in § 12.3.21 placing it on the desired or- 
bit, and switching on its self-gravity. We note that after the 
satellite models are grown to full mass and have become self- 
gravitating, the tidal field of the primary disk galaxy acts to 
truncate their outer parts establishing new tidal radii that are 
smaller than the nominal values in Table [T] However, this 
difference is less than 10% in all cases suggesting that the 
primary disk galaxy model and the host CDM halo Gi have 
similar densities at ri^f. 

To limit computational cost, subhalos were removed from 
the controlled simulations once they reached their maximum 
distances from the disk after crossing, so satellites were not 
permitted to begin a second passage. Uniformly, these dis- 
tances were only slightly smaller compared to the starting 
radii of the orbits 50 kpc). This is because dynamical fric- 
tion is very efficient for these massive satellites causing their 
orbit to decay very rapidly. In all cases, satellites lost > 80% 
of their mass after completing their first passage. This justifies 
our decision to neglect the dynamically insignificant subse- 
quent crossing events. We emphasize that only the self-bound 
core of each satellite was extracted from the simulations at the 
subsequent apocenter and not the unbound material, a signif- 
icant fraction of which remained in the vicinity of the disk. 
Removing the latter would result in potential fluctuations that 
could throw the disk out of its relaxed state, and thus interfere 
with the interpretation of the results. 

We present results for the disk structure after allowing the 
disk to relax from the previous interaction. Consequently, our 
findings are relevant to systems that exhibit no obvious, on- 
going encounters. Due to the complexity of the interaction, 
we determined this timescale empirically by monitoring basic 
properties of the disk structure (e.g., surface density, veloc- 
ity dispersion, thickness) as a function of time. When these 
quantities stopped evolving significantly within radii of inter- 
est (changes of the order of 5- 10% in disk properties were 
considered acceptable), the encounter was deemed complete. 
Typical "settling" timescales were found to be in the range 
^ 500-600 Myr The limiting radius was chosen to be equal 
to 7/?d which contains ^ 99% of the initial mass of the disk. 
It is necessary to adopt some limiting radius because dynam- 
ical times grow with radius and the outer disk regions con- 
tinue to show signs of evolution well after each satellite pas- 
sage. These transient structures can be readily identified in 
the edge-on views of the disk presented in the next section. 

When time intervals between subhalo passages were larger 
than the timescale needed for the disk to relax after the previ- 
ous interaction, we introduced the next satellite immediately 
after the disk had settled from the previous encounter. This 
eliminates the computational expense of simulating the disk 
during the quiet intervals between interactions. 

We compute all disk properties and show all visualizations 
of the disk morphology after rotating the disk to the coor- 
dinate frame defined by the three principal axes of the total 
disk inertia tensor and centering it to its center of mass. The 
motivation behind performing these actions is twofold. First, 
infalling satellites tilt the disk through the transfer of angular 
momentum (see § 13.61 ). Introducing the new coordinate frame 
is important because in the original coordinate frame, rota- 
tion of stars in a tilted disk would appear as vertical motion 
interfering with the interpretation of the results. Of course, 
the tilting angles of the inner and outer parts of the disk may 
differ from each other and from the tilting angle of the disk 
determined by the principal axes of the total disk inertia ten- 
sor. However, these differences are typically small (< 2°; see 



Figure [8]i and we have verified that they only have a minor 
influence on measured disk properties. In what follows, we 
quote all results relative to the coordinate frame defined by 
the three principal axes of the total disk inertia tensor. Sec- 
ond, the masses of the simulated satellites are a substantial 
fraction of the disk mass. As a result, the disk center of mass 
may drift from its initial position at the origin of the coordi- 
nate frame due to the encounters with the infalling subhalos. 

3. RESULTS: DYNAMICAL SIGNATURES OF HIERARCHICAL 
SATELLITE ACCRETION 

In this section we examine the response of the thin, disk 
galaxy model Dl to interactions with cosmological subhalos 
S1-S6 of host halo Gi, which are designed to mimic a typical 
central accretion history for a Galaxy-sized CDM halo over 
the past ^ 8 Gyr The "final" disk discussed in the next sec- 
tions has experienced the S1-S6 encounters and was further 
evolved in isolation for ^ 4.3 Gyr after the last interaction, so 
that the disk evolution is followed from z = 1 to z = 0. We re- 
iterate that though our numerical experiments employ a best- 
fit model for the present-day structure of the MW, our aim 
is to study the general, dynamical evolution of a disk galaxy 
subject to a ACDM-motivated satellite accretion history since 
1. 

The focus of this study is exclusively on the evolution of the 
disk material, so we do not consider the bulge component in 
any of the analysis presented below. Lastly, while our simula- 
tion program mimics the accretion events in halo Gi , the simi- 
larity of subhalo populations in all four Galaxy-sized host ha- 
los we analyzed suggests that the results presented next should 
be regarded as fairly general. We discuss in turn disk global 
morphology, thickening, velocity structure, surface density, 
lopsidedness, and tilting. 

3.1. Global Disk Morphology 

Figure [5] illustrates the global response of the disk to the 
infalling subhalos. The encounter with the first substructure 
(SI) generates a conspicuous warp beyond ^12 kpc. The im- 
pact of the second most massive satellite (S2) has a dramatic 
effect on the global disk structure. The entire disk visually 
becomes considerably thicker compared to the initial model 
after this accretion event. In addition, this interaction excites 
extended ring-like features in the outskirts of the disk and a 
moderately strong bar (Paper I), both of which indicate that 
the axisymmetry of the initial disk is destroyed as a result of 
the first two accretion events. We stress that the bar is in- 
duced in response to the subhalo passages, not by amplified 
noise. It has a semi-major axis varying between ^3-4 kpc, 
within the range of values inferred for the bar in the MW (e.g., 
iBissantz & Gerhardll2002h . Throughout its evolution the bar 
fails to form a boxy bulge and remains rather thin. We note 
that we evolved the disk galaxy in isolation for ^ 5 Gyr af- 
ter the encounter with SI and verified that the first satellite 
alone is not capable of exciting a bar on these timescales. 
Fourier decomposition of the final disk also reveals that the 
radial variation of the amplitude of the observed spiral struc- 
ture is similar to that of normal spiral galaxies in the near 
infrai-ed (e.g., G rosbol & Patsisll 19981) . 

Non-axisymmetric structures such as the bar and the rings 
drive continued evolution b y redistributing mass a nd angular 
momentum in the disk (e.g.. lDebattista et ani2006l) . Thus, in- 
falling satellites affect galactic disks not only directly by im- 
pulsively shocking the orbits of disk stars, but also indirectly 
by exciting global instabilities. In fact, angular-momentum 
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Fig. 3. — Density maps illustrating the global morphological evolution of the disk in galaxy model Dl subject to a typical ACDM-motivated accretion history 
expected for a Galaxy-sized dark matter halo since z ~ 1 . Particles are color-coded on a logarithmic scale with brighter colors indicating regions of higher stellar 
density. Local density is calculated using an SPH smoothing kernel of 32 particles. The face-on (bottom panels) and edge-on (upper panels) distributions of disk 
stars are shown in each frame and to aid comparison the first panel also includes the edge-on view of the initial disk. The labels for individual satellite passages 
from S 1 to S6 and the time coiTesponding to each snapshot are indicated in the upper left-hand and lower right-hand corners of each bottom panel. Results are 
presented after centering the disk to its center of mass and rotating it to a new coordinate frame defined by the three principal axes of the total disk inertia tensor. 
The first satellite passage generates a conspicuous warp while the second encounter, which involves the most massive subhalo S2, causes substantial thickening 
of the disk and excites a moderately strong bar and extended ring-like features in the outskirts of the disk. Accretion histories of the kind expected in ACDM 
models play a substantial role in setting the global structure of galactic disks and driving their morphological evolution. 



transport from one part of the disk to another causes the disk 
to expand radially during the encounters (see § 13. 4b . However, 
we expect this expansion to be significantly less pronounced 
compared to that in the vertical direction as rotational energy 
dominates over random motions in the plane of the disk. 

Apart from changes associated with the structural details of 
the bar and the location and extent of the ring-like features, the 
disk structure does not evolve appreciably during subsequent 
satellite passages (S3-S6). In broad terms, the global morpho- 
logical evolution of the galactic disk is driven by the interac- 



tion with the most massive subhalo of the accretion history. 
We will return to this point repeatedly in subsequent discus- 
sions. Lastly, the face-on projections of the disk show that the 
tidally induced bar and other structures are non-transient fea- 
tures that survive long after the initial perturbations. Indeed, 
we confirmed their presence in the final disk some 4.3 Gyr 
after the last encounter (see also Figure |6] in Paper I). 
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Fig. 4. — Disk thickening. Tliickness profiles, z(i?), of the disk in galaxy model Dl viewed edge-on. Thicknesses and radii are normalized to the scale height, 
Zrf, and radial scale length, Rj, of the initial disk. Left: Thickness profiles for the initial (thin lines) and final disk (tliick lines). Lines of intermediate thickness 
show the fractional increases in disk thickness defined as [median(|z|)f-median(|z|)i]/median(|z|)j, where median(|z|)f and median(|z|)i denote the thickness of 
the final and initial disk, respectively. Dot-dashed, solid, and dotted lines correspond to thicknesses measured by the median of the absolute value, median(|z|), 
the mean of the absolute value, < |z| >, and the dispersion, < >'''^, of disk particle height above the midplane, respectively. The initial disk is constructed 
with a constant scale height, explaining why the corresponding curves are flat. The arrow indicates the location of the solar radius, Rq . The initial thin disk 
thickens considerably at all radii as a result of the encounters with CDM substructure and a conspicuous flare is evident in the final disk beyond R > 4Sj. Right: 
Evolution of disk thickness profiles. Different fines show results for individual sateflite passages from Si to S6 and thicknesses are measured as median(|z|). The 
vertical dotted line indicates R0 and various symbols correspond to the pericenters of the infalling subhalos. The thick solid line corresponds to the initial disk 
evolved in isolation for a timescale equal to that of the combined satellite passages. The first two satellite passages thicken the disk considerably at all radii and 
cause substantial flaring in the outskirts. The combined effect of the remaining subhalos (S3-S6) is much less dramatic, indicating that the second, most massive 
accretion event is responsible for setting the scale height of the final disk. 



3.2. Disk Thickening 

Among the most intriguing dynamical effects of the subhalo 
impacts is the pronounced increase in disk thickness, which is 
evident in the edge-on views of the disk in Figure [3] A more 
quantitative measure of the thickening of the disk in response 
to the satellite accretion events is presented in Figure H) This 
figure shows thickness profiles of the disk, z{R), viewed edge- 
on. 

The thickening of the disk can be formally described by 
the increase of its vertical scale height. In the present study 
we quantify disk thickness at any given radius by the median 
of the absolute value of disk particle height above the disk 
midplane, median(|z|). The motivation behind this choice 
is threefold. First, it permits direct comparison with ear- 
lier work because similar estimators of disk thickness have 
been employed by previous authors (e. g., Ouinn et al .lll993l: 
Walker et al.lll996t Ivelazquez & Whit3Tl999l: iGauthier et alj 
2006^ . Second, the gravitational interaction between satellites 
and disks is capable of heating individual disk stars at quite 
large radii above the plane of the disk. The choice of the me- 
dian value of the particle height mitigates the influence of dis- 
tant outlier particles and thus produces conservative estimates 
for the increase of disk thickness. Third, disk scale heights 
must be formally derived by means of fitting the particle dis- 
tribution to an appropriate functional form (e.g., exponential 
or a sech^ law). In contrast, quantities such as median(|z|) do 
not require a fit to any functional form and so they are unam- 
biguous and require no assumptions for their interpretation. 



The left panel of Figure |4] shows thickness profiles for the 
initial and final disk together with the fractional increase in 
thickness caused by the infalling satellites. For convenience, 
in Fig. m we compare the median(|z|) to two additional es- 
timators of disk thickness used extensively in the literature, 
namely the mean of the absolute value, < |z| >, and the dis- 
persion, < >'''^, of disk particle height above the mid- 
plane. While all different estimators yield nearly identical 
fractional increase in disk thickness, median(|z|) leads to a 
smaller amount of disk thickening at large radii {R > JiRd)- 
This is true for all of our analyses and so we address only 
median(|z|) in the remaining of the paper. 

We note that the initial disk exhibits a small departure from 
a pure exponential profile near its center. This feature devel- 
ops when constructing the composite DF and potential of the 
multicomponent galaxy model. As a result, the stellar disk 
is somewhat thinner in the middle though the resulting den- 
sity is fully self-consistent within the context of the epicyclic 
and third integral approximations that are built into the disk 
DF (Widrow & Dubinski 2005). This deviation is insignifi- 
cant for our purposes because the influence of this effect on 
disk thickness is quite moderate, and is confined to a fraction 
of a scale length R^. The main goal of the present study is 
to assess the global dynamical response of a galactic disk to 
cosmologically-motivated satellite accretion events. 

The left panel of Fig. |4] demonstrates that the initial, thin 
disk thickens considerably at all radii as a result of the inter- 
actions with subhalos S1-S6. More specifically, disk thick- 
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ness near the solar radius has increased in excess of a factor 
of 2. We emphasize that the thickening does not occur uni- 
formly as a function of radius suggesting that the inner and 
outer disk regions respond differently to the accretion events. 
The outer disk is much more susceptible to damage by the ac- 
creting satellites. Indeed, atR = Rd the thickness increases by 
^ 50% compared to a factor of 3 rise at /? ^ ARj. The larger 
binding energy of the inner exponential disk and the presence 
of a massive, central bulge (M/, ~ O.SMdisk) which deepens 
the central potential are likely responsible for the robustness 
of the inner disk. Given the fact that the infalling subhalos 
are spatially extended and the self-gravity of the disk grows 
weaker as a function of distance from the center, it is not un- 
expected that the scale height of the disk should increase with 
radius. Indeed, by making the simplest assumption that the 
accreting satellites deposit their orbital energy evenly in the 
disk, it can be shown that the disk scale height increases as 
Az(7?) cx Tj~^{R), where is the disk surface density (Pa- 

per I). It is also worth emphasizing that a significant part of 
the evolution seen in the inner parts of the disk {R < Rd), in- 
cluding the peak in the thickness profiles at /? ~ 0.5Rd, can be 
attributable to the tidally induced bar. 

The right panel of Figure |4] presents the evolution of disk 
thickness caused by individual satellite passages S1-S6. After 
the SI encounter, the thickness of the outer disk (R > 3.5Rd) 
rises considerably leading to a a distinct flare. Specifically, 
at R = ARd the thickness increases by ^ 60% compared to a 
^ 20% increases at Rq. In contrast, the inner disk (/? < 1 -SRd) 
appears unaffected by the accretion event despite the sig- 
nificant initial mass (Msub — 0.3Mdisk) and small pericenter 

(rpen^l.2/?rf) of SI. 

The second accretion event involves the most massive satel- 
lite S2 (Msub — 0.6Mdisk) and is responsible for generating a 
much more pronounced flare in the outer parts as well as caus- 
ing the entire disk to thicken. However, the inner disk regions 
still exhibit substantial resilience to the encounter Indeed, 
disk thickness at R = \ .5Rd rises only by 40% compared 
to a factor of 3.6 increase at R = ARd- The corresponding 
increase at Rq is a factor of ~ 1 .9. 

Though the masses of satellites SI and S2 differ by less 
than a factor of 2, subhalo S2 constitutes a far more effi- 
cient perturber compared to satellite S 1 . Taking into account 
that both subhalos have similar internal properties (Table [TJ 
and are on approximately the same polar orbit, this suggests 
that disk thickening is very sensitive to perturber mass. In 
fact, us ing A^-body simulations of satellite-disk encounters 
iHayashi & Ch iba (2006) found that the increase of the disk 
vertical scale height,Azf/, is proportional to the square of the 
subhalo mass, Azd/Rd M^^^^. Our results suggest that such 
correlation may be even stronger. 

The remaining satellite passages have a markedly less dra- 
matic effect in perturbing the vertical structure of the disk. 
This fact confirms that subhalo interactions with already 
thickened disks induce much smaller relati ve changes in disk 
thickness compared to the initial passages (lOuinn et al.ll 19931; 
Paper I). It also suggests that thicker disks may exhibit en- 
hanced resilience to e ncounters with substructure. We return 
to this point in § 14.11 The combined effect of subhalos S3- 
S6 acts to increases the disk thickness at intermediate (Rq) 
and large [R = 4Rd) radii by only ^ 10% compared to that af- 
ter passage S2. The larger differences observed at small radii 
(/? < 1 .5Rd) are attributable to the bar, which grows stronger 
as a function of time. Finally, we remark that the thickness of 



the same disk galaxy evolved in isolation for a timescale equal 
to that of all satellite passages grows by only ^ 10% indicat- 
ing the excellent quality of the initial conditions and adequate 
resolution of the simulations. 

3.3. Disk Velocity Structure and Heating 

The dynamical response of a thin, galactic disk to infalling 
satellites manifests itself in a variety of ways. In this sec- 
tion we address the degree to which the disk velocity struc- 
ture evolves as a result of the simulated accretion history. 
We characterize the disk kinematical response to accretion 
events through the evolution of the disk velocity ellipsoid 
(ctr, (70 , 0";), where ctr, ct^, and ct; correspond to the radial, az- 
imuthal, and vertical velocity dispersions, respectively. In the 
isothermal sheet approximation, we expect that <t^ oc Zd, and 
so the evolution of the disk vertical velocity dispersion should 
be less pronounced than the evolution of disk thickness. The 
relevant analysis is presented in Figure |5] 

This figure reveals that the disk velocity ellipsoid grows sig- 
nificantly in all three directions as a result of the gravitational 
interactions with the infalling satellites. The velocity struc- 
ture of the disk evolved in isolation for an equivalent period 
shows virtually no evolution in the vertical direction. Some 
minimal evolution in the radial and azimuthal components of 
the ellipsoid is observed because spiral structure arising from 
the swing amplification of discreteness noise leads to angular- 
momentum transport in the plane of the disk. Nevertheless, 
this increase is negligible compared to the net effect caused 
by the satellite encounters. 

As with disk thickening, the velocity structure is influenced 
most dramatically by S2. Both the direct deposition of en- 
ergy by the subhalos and global instabilities such as the bar 
and spiral structure which are excited during these encounters 
are responsible for the observed evolution. The latter phe- 
nomena act predominantly in the plane of the disk and cause 
the planar components of the disk velocity ellipsoid to evolve 
significantly. 

The disk velocity ellipsoid at the solar radius Rq in- 
creases from (crR,cr0,crj;)= (3 1,24, 17) kms"' to ((TR,a^,a,) ~ 
(6 1, 49, 31) km s"'. For reference, the velocity ellipsoid 
of the thick disk of the MW at Rq is estimated to 
be - (46,50,35) km s"' by IChiba & BeersI (l2000l) and ~ 
(63,39,39) kms"' bv lSoubiran et al.l (l2003h . However, there 
are at least two reasons to exercise caution interpreting the 
details of these results. First, the final dispersions are calcu- 
lated considering all stars, rather than only those that belong 
to the thick disk component (see Paper I). Second, we have ne- 
glected any stellar components in the infalling satellites and 
thus we are unable to address their contribution to the final 
disk. This is relevant in light of theoretical (e.g., Abadi et aLl 
2003) as well as observational studies ( Yoachim & Dalcanto^ 
2006) suggesting that stars formed in accreted satellite galax- 
ies play a significant role in the formation of thick disks. 

Figure |5] shows that subhalo impacts heat the galactic disk 
in a non-uniform way. At small radii {R < Rd), the planar 
velocity dispersions exhibit a steep increase which is again 
primarily a consequence of the bar. The central values of ctr 
and CT0 grow by a factor of ^ 2.2, while <t~ increases by a 
factor of ~ 1 .6. As with disk thickness, the effect of the ac- 
cretion events on the velocity ellipsoid grows stronger with 
distance for R > Rd- While the initial disk is designed so that 
its velocity ellipsoid decreases monotonically as a function 
of projected radius, after the satellite passages all dispersions 



12 



KAZANTZIDIS ET AL. 



I I I I I I I I I I I I I I I I I I I I I I I I I I I I 

initial 
SI 




1 1 1 1 1 1 1 1 1 1 1 1 


1 1 1 1 1 1 1 1 1 1 1 1 1 1 






/ 


^ A. 

\ 


1 1 1 1 1 1 1 1 1 1 1 1 1 1 


1 1 1 1 1 1 1 1 1 1 1 1 1 1 



o 

o 



b 0.5 - 



I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I 




Fig. 5. — Disk heating. Velocity dispersion profiles of the disk in galaxy model Dl as a function of projected radius in units of the radial scale length of the 
initial disk, R^. Counterclockwise from the upper left, the panels display the evolution of the radial, or, azimuthal, cr^, and vertical, u-^, velocity dispersions. 
All profiles are normalized to the total circular velocity of galaxy model Dl at the solar radius, Vc(^q) = 234.1 kms"'. Different lines con'espond to individual 
satellite passages from SI to S6 and line types are the same as in the right panel of Figure |4] Thick lines present the kinematical properties of the primary 
disk galaxy evolved in isolation for a timescale equal to that of the combined satellite passages. In the upper right panel the velocity dispersion ratios cr-/a^ 
(solid lines) and u^jdR (dot-dashed lines) are shown for the initial (thin lines) and final (thick lines) disk. In all panels arrows indicate the location of the solar 
radius, Rq. Bombardment by CDM substructure heats the thin galactic disk considerably in all three directions and causes its velocity ellipsoid to become more 
anisotropic. 



become nearly constant outside ^ 3/?,/. Such flat velocity dis- 
persion profiles are in very good agreement with results from 
recent kinematic studies o f planetary nebulae in the extreme 
outskirts of spiral galaxies dHerrmann et al.ll2009h . 

In addition, though the initial dispersion profiles are 
smooth, at radii R > R^, the planar components or and 
of the velocity ellipsoid exhibit wave-like features which are 
associated with the stellar rings in the disk plane seen in Fig- 
ure|3] These structures result from the satellite encounters and 
do not dissipate even after several Gyr of evolution subsequent 
to the last accretion event. 

The upper-right panel of Figure |5] shows that the ratio of 



(T, to decreases at all radii as a result of the subhalo im- 
pacts. This finding, in conjunction with the fact that ct^/o-r 
appears not to be affected in any noteworthy way, indicates 
that both planar components of the disk velocity ellipsoid re- 
spond more strongly to the accretion events compared to the 
vertical component. In other words, infalling satellites cause 
the disk velocity ellipsoid to become more anisotropic. The 
larger increase of ctr and relative to a, is due to the pres- 
ence of the bar and the action of spiral structure that transports 
angular momentum through the disk. It is worth noting that 
though the resultant heating is substantially larger in the plane 
the disk spreads primarily in the vertical direction. This is be- 
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Fig. 6. — Disk antitruncation. Evolution of tlie surface density profiles, 
S(/?), of the disk in galaxy model D 1 viewed face-on as a function of pro- 
jected radius, R, in units of the radial scale length of the initial disk, Rj. 
Different lines show results for both the disk initially (solid line) and the disk 
after each subhalo impact as indicated in the upper right-hand corner. All 
profiles are normalized to the surface density of the initial disk at the origin. 
So, and the vertical dotted line indicates the location of the solar radius, Rq. 
Encounters with CDM substructure can generate surface density excesses in 
galactic disks similar to those seen in the light profiles of observed antitrun- 
cated disk galaxies. 

cause rotational energy dominates over random motions in the 
plane of the disk. 

3.4. Disk Surface Density and Antitruncation 

In is interesting to investigate how the simulated subhalo 
accretion history would influence the disk surface density dis- 
tribution. Figure |6] shows the evolution of t\\t face-on surface 
density profile of the disk as a function of satellite passages 
S1-S6. By construction, the surface density profile of the ini- 
tial disk follows an exponential distribution in cylindrical ra- 
dius R. 

As with both disk thickness and velocity structure, the disk 
surface density is minimally affected by the encounter with 
the first satellite S 1 . The response of the surface density distri- 
bution to subsequent accretion events is notable and in agree- 
ment with the results of the previous subsections it is driven 
by the most massive subhalo S2. 

At small radii {R < Rj), the surface density profile steep- 
ens considerably, a feature that is again due to the bar and 
the transport of disk material inwards by the spiral patterns 
and ring-like features. The central disk surface brightness in- 
creases by A/i = 1 mag arcsec"^ during the course of the inter- 
actions with subhalos S1-S6. At intermediate radii {R^ ^R ^ 
5Rd), the surface density profiles are very similar to that of the 
initial disk. Additionally, these profiles are not smooth but 
display strong wave-like features which are associated with 
the rings of disk material (Figure[3]l. 

More interestingly, the disk spreads radially, and beyond 
R > 5Rd, there is a clear excess of surface density compared 
to that of the initial disk. The subhalo encounters modify 
the face-on structure of a single-component, exponential disk 
to producing a system with two distinct components. Be- 



yond the very inner disk regions where the bar dominates 
and within approximately 5 scale lengths of the initial disk 
(5Rd ^ 14 kpc), the final surface density is exponential with 
a scale length roughly equal to that of the initial disk. In- 
deed, fitting an exponential profile to the mass surface den- 
sity of the final disk at Rj < R < 5Rd yields a best-fit scale 
length of Rd, inner ~ 3.05 kpc, or ~ 8% larger than the initial 
scale length. At larger radii the disk exhibits a flatter surface 
density profile, which is again crudely exponential. Fitting 
another exponential profile to the outer final disk (R > 5Rd) 
yields a scale length of Rd,outer ~ 5 kpc ^ l.^Rd ^ i.6Rd.imKi- 
This excess surface density at large radii relative to the ex- 
ponential profile of the inner disk is interesting in the con- 
text of the so-called "antitruncated" stellar disks whose sur- 
face brightness profiles display similar excesses at roughly 
4 — 6 disk scale lengths from galaxy centers (e.g., Erwin et al 



2005tlPohlen & Truiillol2006[lPohlen et al. 120071; Erwin et al 



2008D. 

Our simulations lack the physical effect of star formation, 
so this disk antitruncation arises from old stars that resided 
in the initial thin disk and migrated outward. The driver be- 
hind the expansion of the disk and the excess surface density 
is redistribution of disk angular momentum (and hence disk 
mass) caused both directly and indirectly by the subhalo im- 
pacts. Infalling satellites directly deposit kinetic energy and 
angular momentum into the disk during the gravitational in- 
teractions. Additionally, non-axisymmetric features such as 
the bar and spiral structure transport angular momentum and 
stellar mass to large radii. In response to the net transfer of an- 
gular momentum between its inner and outer regions, the disk 
expands. Indeed, we find that the magnitude of the total angu- 
lar momentum at 7? > 5Rd grows by ^ 80% over the course of 
the simulated satellite accretion history. The transport of disk 
material outwards in radius also leads to the excess surface 
density observed at large projected radii. 

3.5. Disk Lopsidedness 

The face-on views of the simulated disk in Figure |3] illus- 
trate that encounters with infalling satellites are responsible 
for generating large-scale, long-lived asymmetries. An 
intriguing phenomenon observed in many disk galaxies is the 
existence of significant lopsided asymmetries in their neutral 
hydrogen and/or stellar mass distributions (e.g., Baldwin et aH 
1980; Richter & Sancisi 1994; Rix& Zaritskv 19951 
Zaritsky & Rix 1997; Matthews et al. 1998; Havnes et aQ 
179981: IBournaud et al . 2005; ReichardetalJL2008).' System- 
atic attempts to quantify the frequency of such asymmetries 
in stellar disks have revealed the presence of lopsidedness 
in a substantial fraction of disk galaxies. The fraction of 
lopsided disks varies between ~ 20% and 30% for strongly 
lopsided systems and may exceed 50% for moderately lop- 
sided disks (e.g., Rix & Zaritskv 1995; Zaritskv & Rix 19971 
iRudnick & R ix 1998; Bournaud et al. 2005). It is therefore 
interesting to quantify the degree to which the simulated disk 
exhibits deviations from axisymmetry. 

A standard characterization of asymmetries in galaxies is 
via Fourier decomposition. We express the surface density of 
the stellar disk as a Fourier series 



/m[0-0„,(«)] 



(3) 



m=l 



where A,„(/?) = a,n{R) / aQ{R) and 4>i„{R) denote the normalized 
strength and the phase of the Fourier component m, and quan- 
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tify lopsidedness by the ratio of the amplitudes of the m = I 
to m = (azimuthally-average d surface density) Fou rier coef- 
ficients, A i (7?) ee a i(/?)/flo(/?) (" Rix & Zar itskv'l995h. We per- 
formed the Fourier decomposition after centering the system 
to the peak of the density distribution (Rix & Zaritskv" 19951) . 
whereas using centroids ( Debattista & Sellwood 200Q) gave 
similar results. 

Figure |7] shows the variation of the lopsidedness parame- 
ter A i as a function of projected radius from the center of the 
disk. Results are presented for both the initial and final dis- 
tribution of disk stars and for four characteristic timescales in 
the simulated accretion history of halo Gi. The initial disk 
is constructed to be axisymmetric explaining why the corre- 
sponding curve is flat. Fig. [T] clearly demonstrates that en- 
counters with CDM substructures trigger lopsidedness in the 
galactic disk. Furthermore, the induced lopsidedness is not 
constant as a function of radius nor are changes monotonic. 
Different regions of the disk may thus become lopsided to 
different degrees by the infalling satellites. 

During the simulated accretion history, typical values of 
Ai span the range 0.1 < Aj < 0.2 (l.5Rd < R < 4.5Rd), but 
more prominent lopsided asymmetries reaching amplitudes of 
~ 0.3 is observed at larger radii (R > 5Rd)- These values are 
consistent with observational estimate s of lopsidedness from 
various samples of stellar disks (e .g., 'Rix & Zaritskv"1995t 
Zaritskv & RIxI|1997I: iRudnick & R ix 1998; Bournaud et al. 
20051) . The outer disk regions are more susceptible to the 
tidal perturbations that generate these asymmetries, and be- 
cause they are characterized by longer dynamical times, sig- 
nificant lopsidedness persists there. Though the bombard- 
ment by CDM substructure has ceased by 3.5 Gyr, lopsid- 
edness of moderate amplitude is still imprinted in the stellar 
disk ^ 2 Gyr after the last accretion event. The final disk 
which corresponds to ~ 4.3 Gyr after the last satellite im- 
pact exhibits virtually no signs of lopsided asymmetries at 
R < 4.5Rd and very weak lopsidedness at larger radii. The 
lifetimes of the reported asymmetries (~ 1 Gyr) are in agree- 
ment with e stimates from phase mixing and winding argu- 
ments (e.g.. [Baldwin et al.|[T980t iRix & Zari tskv" 1995) and 
direct analysis of A^-body s i mulations (e.g ., Zaritskv & Rix 
[19971 iBoumaud et aP ^2m5; MapeUieialJ[2008). A more 
thorough discussion regarding the longevity of lopsidedness 
excited by halo substructure is beyond the scope of the present 
paper and is deferred to future work. 

The main implication of the results reported in Figure |7] 
is that the continuous accretion of satellites over a galaxy's 
lifetime constitutes a significant source of external pertur- 
bations that excite as well as maintain lopsidedness in stel- 
lar disks at observed levels. This might be important as 
many observational studies find no correlation between the 
presence of nearby companions and lopsidedness in disk 
galaxies (e.g., Zaritskv & Rix.l997:.Wilcots & Prescott.2004; 
iBournaud et al.ll20051) ". 

3.6. Disk Tilting 

As discussed earlier, angular momentum exchange between 
the satellites and the disk is expected to tilt the disk relative 
to its initial orientation. The left panel of Figure [8] presents 
edge-on density maps of the initial and final disk. The final 
disk is shown relative to both the initial coordinate frame and 
the frame defined by the principal axes of the total disk inertia 
tensor 

The panel reveals that the simulated subhalo impacts cause 
a significant amount of disk tilting. For this particular accre- 
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Fig. 7. — Disk lopsidedness. Evolution of the Ai parameter, the normalized 
ampHtude of the m = i Fourier component of the disk surface density, of the 
disk in galaxy model Dl as a function of projected radius in units of the 
radial scale length of the initial disk, Rj. Results are presented for both the 
initial {open circles) and final distribution of disk stars (filled circles), and 
for four characteristic times in the simulated accretion history of host halo 
Gi (after each of the SI, S3, and S6 encounters, and ~ 2 Gyr subsequent 
to the last accretion event). ACDM-motivated satellite accretion histories 
are responsible for triggering as well as maintaining for a significant fraction 
of the cosmic time lopsidedness in stellar disks at levels similar to those in 
observed galaxies. 

tion history, the angle by which the disk is tilted from its orig- 
inal plane is ^ 1 1°, implying a significant transfer of angu- 
lar momentum from the infalling satellites to the host galactic 
disk. This indicates that not all of the orbital energy associ- 
ated with the vertical motion of the subhalos is converted into 
random motions of disk stars causing disk thickening. When 
viewed in the original coordinate frame, the final disk appears 
visually much thicker, which is simply a consequence of the 
fact that the disk is substantially tilted by the orbiting satel- 
lites. Employing the tilted coordinate frame for all computa- 
tions of disk properties is required as, in the original coordi- 
nate frame, quantities such as the scale height, z{R), and the 
vertical kinetic energy, T, will be artificially inflated. Before 
continuing, we note that in the absence of satellites there is 
some exchange of angular momentum between the disk and 
the dark matter halo of the primary galaxy. This exchange re- 
sults in tilts of < 1 ° when the galaxy is evolved in isolation for 
a timescale equal to that of the combined subhalo passages. 

The right panel of Figure [8] shows a scatter plot of the evo- 
lution of the angular position of disk pole, defined by the az- 
imuthal and zenith angles (cj), 9), relative to initial disk pole 
as a function of satellite passages S1-S6. These angles are 
computed considering all disk particles within 15 kpc from 
the disk center, but we note that the results do not depend 
sensitively upon this cutoff. The combined action of the first 
two satellites S 1 and S2 serves to tilt the galactic disk by a 
very small amount 6 < 1.5°. Both subhalos are on nearly po- 
lar orbits (Table [T]) and thus transfer of angular momentum 
to the disk is expected to be minimal. In this case, virtu- 
ally all of the orbital energy associated with the z motion of 
the infalling satellites is converted into random vertical stellar 



T — I — I — I — I — I — I — I — I — I — I — I — I — I — |~i — I — I — I — I — I — I — I — r 
O initial disk 



tsi = 0.6 Gyr 
- ■ ts6=3.4 Gyr 




CDM Substructure and Galactic Disks II: 



15 




■VAMBBMMMH H^JL^^ KMMMa KMMWB KMSMflMMMMB KM^£ 



CD 6 - 



I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I 




I I I 





90 180 270 360 

0[1 



Fig. 8. — Disk tilting. Left: Density maps of the disk in galaxy model Dl viewed edge-on. The upper panel shows the initial disk, whereas the middle and 
bottom panels depict the final disk in the initial coordinate frame and the frame defined by the three principal axes of its total inertia tensor, respectively. The 
zenith angle 8 between the two coordinate frames is indicated in the middle panel. Right: A scatter plot of the evolution of the angular position of disk pole 
relative to initial disk pole, where c/) denotes the azimuthal angle. Different symbols correspond to individual satellite passages from SI to S6. Filled circles show 
results for the initial disk evolved in isolation for a timescale equal to that of the combined satellite passages. The inset presents the evolution of angle S as a 
function of radius from the center of the disk in units of the disk truncation radius, Sout = 30 kpc. Interactions with infalling satellites of the kind expected in 
ACDM models can drive substantial tilting in galactic disks. 



motions causing thickening of the disk. The next four inter- 
actions involve subhalos on both prograde (S3,S5) and ret- 
rograde (S4,S6) orbits, so angular momentum exchange can 
take place. While these impacts result in 6' ^ 11°, they keep 
between (0-35)°, so they tilt the disk in almost the oppo- 
site direction compared to passages SI and S2. It is interest- 
ing that these four accretion events continue to tilt the disk 
in this small range of 0. Lastly, the disk response to subhalo 
passages S3-S6 suggests that, in addition to the mass, the or- 
bital orientation of the infalling satellite is crucial in determin- 
ing the amount of disk tilting. Retrograde encounters appear 
to be associated with more pronounced tilting compared to 
thei r pro grade counterparts. We address this issue explicitly 
in §133] 

The inset in the right panel of Figure [8] presents the evo- 
lution of the angle 6 between the original and tilted coordi- 
nate frames as a function of radius from the center of the 
disk. The maximum radius we consider for this calculation 
is r = 0.85/?out = 25.5 kpc, where /?out denotes the disk trun- 
cation radius ( Widrow & D ubinski 2005), which is set by the 
requirement that each radial bin should contain at least 1000 
particles. This is an empirical criterion which ensures a ro- 
bust determination of the inertia tensor in each bin. The inset 
illustrates that the angle 6 between the original and tilted co- 
ordinate frames is not constant as a function of radius and 
changes are not monotonic. Different regions of the disk may 
thus be tilted to different degrees by the infalling satellites. 
Interestingly, the dense inner disk (r/7?out ^ 0.5), which con- 
tains most of the disk mass, r esponds nearly as a rig id body to 
the accretion events (see also lShen & Sellwoodll2006i) . While 
the tilting angles of the inner and outer parts of the disk are 
different from each other, which is indicative of the pres- 
ence of warping, the maximum difference is relatively small 
(A6' < 2°). Furthermore, as inertia increases quadratically 



with distance, the direction of the axes in the tilted coordi- 
nate frame is mainly determined by the outer parts of the 
disk. These facts justify our decision to employ the coordi- 
nate frame defined by the principal axes of the total inertia 
tensor of the tilted disk to compute disk properties rather than 
to adopt the appropriate local coordinate frame for each radius 
within the tilted disk. 

4. SENSITIVITY OF THE DISK DYNAMICAL RESPONSE TO 
MODELING CHOICES 

In this section we explore various factors that could in- 
fluence the dynamical response of galactic disks to infalling 
satellites. We discuss in turn initial disk thickness, the pres- 
ence of a bulge component in the primary disk, the internal 
density distribution of the infalling systems, and the relative 
orientation of disk and satellite angular momenta. 

4.1. Initial Disk Thickness 



The results in § l3.2l and §[33] suggest that subsequent sub- 
halo encounters with already thickened disks produce smaller 
changes in disk thickness and disk velocity ellipsoid com- 
pared to the initial accretion events. In what follows we in- 
vestigate the effect of initial disk thickness on the thicken- 
ing and heating of a galactic disk by infalling satellites in a 
more controlled manner. For this reason we constructed a 
self-consistent disk galaxy identical to the standard thin disk 
model but with a scale height that was larger by a factor of 
2.5, Zd = 1 kpc (model D2). We then repeated the impacts of 
satellites S 1 and S2 onto this thicker disk model. 

Figure|9]shows the evolution of disk thickness as quantified 
by median(|z|), and disk velocity ellipsoid {<jr,(j^,(jS), after 
each subhalo passage. There is little additional thickening in 
model D2 within roughly four scale lengths. Even beyond 
this radius, D2 is much less flared compared to disk model 
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Fig. 9. — The effect of the initial disk thickness on the response of galactic disks to encounters with satellites. Bottom left panel: Thickness profiles as a 
function of projected radius in units of the radial scale length of the initial disk, Rj. Thin lines show the evolution of thickness in the standard disk galaxy model 
Dl {zd = 0.4 kpc) after the interactions with subhalos SI and S2. Thick lines present results for the same impacts onto the much thicker disk galaxy model 
D2 {zd = 1 kpc). Dashed and dotted Hues correspond to satellites SI and S2, respectively, while the solid curves show results for the initial disks. All curves 
are normalized to the initial scale height of each corresponding disk. Upper left panel: Thickness increase, Az = [median(|z|)(inai - median(|z|)i„itiai]. Here 
median(|z|)ini[ii,i and median(|z|)finai denote the thicknesses of the disk initially and after each subhalo impact, respectively. Right: Evolution of the disk velocity 
ellipsoid in galaxy model D2. Bottom to top: a-, rj^ + 100 kms"' , and ctr + 175 kms"'. The Une types are as in the left panel, but the results for galaxy model 
Dl are not shown. All profiles are normalized to the total circular velocity of model D2 at the solar radius, V^(Rq) = 235.6 kms"' . The vertical dotted lines in all 
panels indicate the location of the solar radius, Rq. Thicker disks are less susceptible to damage by infalling subhalos compared to their thin counterparts. 



Dl. More specifically, D2 thickens by less than 50 pc at Rq, 
while the thickness of Dl increases by more than 200 pc at 
the same radius. Likewise, we find that the velocity structure 
of the thicker initial disk is significantly less influenced by the 
perturbers. The velocity ellipsoid of model D2 increases by 
less than ^ 15% at the solar radius due to the interactions with 
subhalos SI and S2. The corresponding increase for galaxy 
model Dl is slightly less than a factor of 2 (Figure |5]l. It is 
also worth noting that model D2 does not form a bar as a 
result of the subhalo impacts and is associated with a much 
less prominent spiral structure compared to the thinner disk 
Dl. These findings confirm that thicker disks exhibit a much 
weaker global dynamical response to accretion events com- 
pared to their thin counterparts. 

4.2. Presence of a Bulge Component 

The results presented in Figure|4]demonstrate that the inner 
disk regions exhibit substantial resilience to encounters with 
subhalos. Apart from the larger binding energy of the inner 
exponential disk, the presence of a massive, central bulge with 
(Mh — 0.3M(jisk) in disk model Dl may be responsible for the 
robustness of the inner disk. To address the effect of a bulge 
in reducing the damage induced in galactic disks by infalling 
satellites, we repeated encounters SI and S2 after adiabati- 
cally evaporating the bulge component from disk model Dl 
(see § 12.3.3b . We refer to the resulting bulgeless disk galaxy 
model as D3. 

Figure [To] shows the evolution of disk thickness after each 
subhalo passage for both galaxy models Dl and D3. Not sur- 
prisingly, this exercise demonstrates that a massive bulge en- 
hances the robustness of galactic disks to accretion events. 
Analysis of the disk velocity ellipsoids supports this con- 
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Fig. 10. — The effect of a bulge component on disk thickening. Thin lines 
show thickness profiles for disk galaxy model Dl after the encounters with 
subhalos SI and S2. Thick lines present results for the same impacts onto 
the bulgeless disk model D3. Thicknesses and radii are normalized to the 
scale height, Zd, and radial scale length, Rd, of the initial disk. Dot-dashed 
and dashed lines correspond to satellites SI and S2, respectively, while the 
solid curve shows results for the initial disk. A massive bulge enhances the 
resilience of galactic disks to satellite impacts by both absorbing part of the 
orbital energy deposited by the infalling subhalos and stabilizing the disks 
against the development of global non-axisymmetric instabilities. 
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elusion. As expected, the relative effect of the bulge grows 
weaker as a function of radius from the center of the disk. 
The combined action of satellites SI and S2 increases the 
thickness of model D3 at Rq by a factor of ^ 2.6, compared 
to a factor of ~ 2 in the case of model Dl. At larger radii 
(R > 3.5Rd), model D3 is characterized by a more distinct 
flare compared to Dl. This is interesting as the bulge is lo- 
cated in the inner regions of the disk and suggests that global 
instabilities which are curtailed by a massive bulge are effi- 
cient at driving significant evolution in the outer disk. In- 
deed, we stress that subhalo passages S 1 and S2 excite a much 
stronger bar in model D3 compared to Dl and that this bar al- 
ready develops in response to the interaction with S 1 . 

4.3. Satellite Internal Density Distribution 

It is worthwhile to examine the importance of the satel- 
lite internal density distribution to the structural changes in- 
duced in galactic disks. This interest is motivated by the po- 
tential differences in subhalo structure in models with mod- 
ified dark matter properties and density fluctuation spe ctra 
(e.g.. iHogan & Dalcantonll2000t lAvila-Reese et al.ll20()lh . or 
the presence of baryonic components in some fraction of the 
accreting systems. To this end, we repeated the encounters 
between disk galaxy model Dl and satellites SI and S2 after 
assigning different density profiles to each infalling subhalo. 

We employed density profiles that differ only in their 
asymptotic inner slopes, 7 [see Eq. ([T]i] and required the total 
bound mass of each object to be fixed to its fiducial value. We 
studied two, relatively extreme options for the satellite density 
distributions with the intent that these bracket the range of po- 
tential subhalo structures. The first corresponds to a "steep" 
profile with an inner power-law index equal to 7 = 2, while 
the second follows a "shallow" profile with a constant density 
core, 7 = 0. The size of the density core was chosen to be 
approximately equal to 10% of the tidal radius of each satel- 
lite, > 2 kpc (Table [U. The steep profile addresses the pos- 
sibility of a subhalo containing a galaxy and which may have 
undergo ne adiabatic contraction in response to the growth of 
baryons (Blum enthal et 31.11198 6*), while the cored profile is 
designed to overestimate any density suppression that may re- 
sult in modified dark matter models. The density, p(r), and 
cumulative mass profiles, M(r), for all initial satellite models 
are presented in the left panel of Figure [TT] 

The right panel of Figure [TT]presents the evolution of disk 
thickness in these alternative models compared with the stan- 
dard cases. Not surprisingly, the amount of disk thickening 
increases with the steepness of the inner satellite density pro- 
files at fixed initial subhalo mass. The combined action of 
cored satellites S 1 and S2 increases the disk thickness at the 
solar radius by ~ 45%, compared to a factor of ^ 2 in the 
case of the standard satellites. The cumulative effect of the 
steep subhalos S 1 and S2 is to increase the disk thickness by 
nearly a factor of 2.5 at Rq. The relative thickness dif- 
ferences driven by the different density distributions becomes 
more pronounced at large radii where the disk self-gravity, 
and thus its restoring force, is weaker. It is worth noting that 
steep satellites S 1 and S2 excite a stronger bar compared to 
their 7=1 counterparts, whereas cored satellites S 1 and S2 
do not induce a bar at all on the timescales of these simula- 
tions. These findings indicate that satellite density structure as 
characterized by the inner slope or concentration is an impor- 
tant factor in determining the amount of damage and global 
dynamical evolution that subhalos impart upon disks. This 
is an important point in light of the fact that most previous 



numerical investigations of satellite-disk interactions did not 
realize their satellite models with the exact density structure 
of cosmological subhalos. 

The primary reason for the dependence of disk dynami- 
cal evolution on satellite structure is that subhalos with shal- 
lower profiles are significantly less tightly bound. As a con- 
sequence, they generate much smaller potential fluctuations 
when they cross the disk and correspondingly cause less thick- 
ening and heating. In addition, cored satellites lose mass more 
efficiently even prior to the actual disk crossing. As a result, 
they penetrate the disk with less bound mass, and thus less en- 
ergy and angular momentum available to be delivered to the 
disk stars. Owing to their lower self-gravity, the outer disk 
regions are expected to be more sensitive to the amount of or- 
bital energy deposited by the infalling satellites. This is also 
confirmed by the findings presented in the right panel of Fig- 
ureE] 

4.4. Orientation of Satellite Orbits 

It is interesting to investigate any correlation between disk 
dynamical response and the orientation of satellite orbits. To 
this end, we repeated satellite passages SI and S2 varying 
the angle 9 between their orbital angular momentum and the 
initial angular momentum of the disk in galaxy model D 1 . We 
considered two cases with identical initial orbital inclinations 
with respect to the plane of the host disk (/ = 45°), but exactly 
opposite directions with respect to its rotation. The first is 
a prograde orbit with 6* = 45° and the second is a retrograde 
orbit with 0= 135°. 

Figure [12] compares the results of these simulations with 
our fiducial experiments in which S 1 and S2 are on nearly po- 
lar orbits with 9 ~ 90° (Table [T]i- Encounters with subhalos 
on prograde orbits are more efficient at thickening the disk 
compared to their retrograde counterparts. This is suggestive 
of a resonant coupling between the satellite and disk angular 
momenta that is suppressed in retrograde encounters. After 
the completion of passage S2, differences in disk thickness 
between prograde and retrograde interactions manifest at all 
radii. Prograde satellites SI and S2 increase the disk thick- 
ness by ^ 20% more compared to the corresponding retro- 
grade models. Moreover, retrograde SI and S2 impacts excite 
a weaker bar than the prograde encounters. 

Interestingly, neither prograde nor retrograde orbits cause 
as much thickening (or heating) as the nearly polar orbits of 
SI and S2 in our fiducial simulations. However, the combined 
action of SI and S2 tilts the disk hy 9 ^ 7.5° in the prograde 
case and 6 ^ 9.5° in the retrograde case. This is to be com- 
pared to 9 < 1.5° in the fiducial case where SI and S2 are on 
nearly polar orbits. This reinforces our interpretation in § 13.61 
of the different degrees of tilting among encounters S1-S6 as 
a manifestation of the different orbital orientations of the in- 
falling satellites relative to the disk. During the encounters, 
part of the orbital energies and angular momenta of the in- 
falling satellites are transferred into the disk. A fraction of 
this energy is converted into random vertical motions of disk 
stars, causing disk thickening, and a fraction is added to the 
disk coherently leading to its tilting. Because polar orbits do 
not transfer angular momentum, the fraction of the satellite 
orbital energy that thermalizes in the disk causing its thicken- 
ing is larger On the other hand, efficient exchange of angu- 
lar momentum results in a large-scale, coherent tilting of the 
disk, and thus in smaller degrees of thickening, an effect that 
is most pronounced in the retrograde orbits. 

We note that since we only consider interactions with satel- 
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Fig. 1 1. — The effect of satellite internal density distribution on disk thickening Left: Spherically-averaged density, p{r), {thin lines) and cumulative mass 
profiles, M(r), {thick lines) for infalling subhalos with different internal structure as a function of radius in units of the tidal radius, rtij. Densities and masses are 
normalized to the enclosed density within the tidal radius, ptid = 3Mbound/47rrj'jj and the bound mass of the system, Mbound {right axis), respectively. Results are 
presented for cosmological satellite S2. Density profiles differ only in the asymptotic inner slope 7 in Eq. (T). The short dashed - long dashed lines show results 
for the standard density profile with 7=1, while the dotted and dot-dashed lines correspond to profiles with an inner power-law index equal to 7 = and 7 = 2, 
respectively. The adopted density disfiibutions have significantly different shapes, but they correspond to exactly the same bound mass. Right: Evolution of disk 
thickness induced by satellites SI and S2 following the density profiles described above. Thin and thick lines display results for the passages of satellites SI and 
S2, respectively, and line types are the same as in the left panel. Solid line shows results for the initial disk model Dl. Thicknesses and radii are normalized to 
the scale height, Zd, and radial scale length, R^, of the initial disk. The vertical dotted line indicates the location of the solar radius, Rq. Infalling subhalos with 
steeper density distributions are more efficient perturbers compared to their cored counterparts. 
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Fig. 12. — The effect of satellite orbital orientation on disk thickening. 
Evolution of disk thickness profiles for disk galaxy model Dl in response to 
encounters with subhalos SI (thin lines) and S2 {thick lines). Thicknesses 
and radii are normalized to the scale height, Zd, and radial scale length, R^, 
of the initial disk. Dashed lines correspond to our fiducial experiments with 
d ~ 90° . Dot-dashed and dotted lines show results for satellites S 1 and S2 on 
aprograde (d = 45°) and a retrograde orbit (9 = 135°), respectively. Solid line 
shows results for the initial disk model D 1 . The vertical dotted line indicates 
the location of the solar radius, Rq. Prograde encounters are more efficient 
at thickening the disk compared to their retrograde counterparts. 



lites that each pass the disk only once, our numerical exper- 
iments cannot address the importance of satellite decay rate 
and mass loss in affecting any of the aforementioned trends. 
This is especially relevant in the case of retrograde and po- 
lar orbits which are known to suffer slower tidal disruption 
and orbital decay compar ed to their prograde counterparts 
(IVelazquez & Whitel[T99l . 

5. DISCUSSION 

We have demonstrated that encounters with halo substruc- 
ture in the context of the ACDM cosmogony imprint a wealth 
of dynamical signatures in thin galactic disks. These signa- 
tures include considerable thickening and heating at all radii, 
surface density excesses resembling those of observed an- 
titruncated disks, lopsidedness at levels similar to those mea- 
sured in observational samples of disk galaxies, and signif- 
icant tilting. In Paper 1, we also showed that the same ac- 
cretion events produce conspicuous flares, bars, low-surface 
brightness ring-like features in the outskirts of the disk, faint 
filamentary structures above the disk plane, and a complex 
vertical structure that is well-described by a superposition of 
thin and thick disk components. All of these findings high- 
light the significant role of encounters with CDM substruc- 
ture in setting the structure of disk galaxies and driving galaxy 
evolution. 

In the present study, we have only utilized a conser- 
vative subset (0.2Mdisk A^sub ^ A^disk, '"peri ^ 20 kpc) of 

a typical accretion history of a Galaxy-sized host halo to 
seed our controlled satellite-disk encounter simulations, and 
have neglected interactions with extremely massive subhalos 

i ^sub ^, ^disk) that could provc ruinous to thin-disk survival 
Purcell etal.ll2008h . In this sense, our results are relevant 
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to systems that have already experienced the most destruc- 
tive events since z = 1 and have re-grown their thin disks 
since. Extending our selection criteria to include subhalos 
with smaller masses (Msub ^ 0.2Mdisk) and/or larger pericen- 
ters (rperi > 20 kpc) will definitely result in many more accre- 
tion events to consider This will perhaps change the detailed 
results for the disk dynamical response subject to a ACDM 
accretion history, but it will not affect any of the qualitative 
conclusions of this study. If anything, the tidal effects of halo 
substructure will be more pronounced compared to that we 
reported here. 

We also emphasize that our numerical experiments were not 
designed to elucidate the effect of subhalo bombardment on 
the structure of specific galaxies such as the MW or M31. 
Rather, the main goal of the present study was to investi- 
gate the most generic dynamical signatures induced in thin 
galactic disks by a typical ACDM-motivated satellite accre- 
tion history. To this end, we utilized only four Galaxy-sized 
dark matter halos to extract representative orbits and subhalo 
merger histories, and the controlled simulations of satellite- 
disk interactions were based on the accretion history of just 
one of these host halos. Nevertheless, all four of these host 
halos showed similar numbers of substantial accretion events 
onto their central regions, and their acc retion histories are typ - 
ical of systems in this mass range (e.g..l Wechsler et al.ll2002h . 
In addition, the bulk of the disk dynamical response illus- 
trated in Figures ID |5] |6] and |7] was driven by the single 
most massive substructure of the simulated accretion history 
(Msuh ^ 0.6Mdisk, /"tid — 20 kpc, and rperi < 10 kpc). We iden- 
tified infalling satellites of this kind in all four of the host 
halo histories studied (Figure [T]!, and such accretion events 
should have been ubi quitous in the histor y of Galaxy-sized 
halos since z ^ 1 (e.g.. lStewart et al]|2008h . All of these facts 
suggest that our simulation set has reasonably captured the 
global dynamical effects of halo substructure on thin galactic 
disks. For detailed comparisons with specific disk galaxies it 
would be required to explore a range of accretion histories and 
orbital distributions of infalling objects using a larger sample 
of halos from cosmological A^-body simulations. 

The analysis presented in Figures|4]and|5]highlight two dis- 
tinctive signatures that infalling satellites imprint in the struc- 
ture and kinematics of the host galactic disk. Figure|4]demon- 
strates that the thickness of the simulated disk increases with 
galactocentric radius, while Figure |5] illustrates that the disk 
velocity dispersion profiles become nearly flat at large radii 
(> 3Rd)- Both features appear to be a natural result of encoun- 
ters with CDM substructure and have important observational 
consequences. The absence of these signatures in a significant 
fraction of disk galaxies would be difficult to reconcile in the 
context of the present study and could be used to potentially 
falsify our proposed model. 

Interestingly, disk-flaring is seen in the MW in both 
the stellar disk ( Lop ez-Corredoira et al.ll2002L iMo manv et al.' 



2006 ) and atomic hydrogen d istribution (e.g., Mer rifield 



1992; iNakanishi & Sofu^ 2003h . Flaring is also observed 
in edg e-on, external galaxies, in their stellar light (e.g., 
Ide Gr£ s & Peletier 1997; Naravan & Jog 2002) and Hj gas 
(e.g..lBrink s & Burton 1984; OlHng 1996; Matthews & WoodI 
120031) . Flat vertical velocity dispersion profiles are also 
reported in the out skirts of disk galaxies. Recently, 
[Herrmann et al.l (12009 ) performed a kinematic study of plan- 
etary nebulae in the extreme outer disks of the nearby, nearly 
face-on spirals M83 and M94. In both these systems, the 
kinematic evidence suggests that; (1) the stellar disks flare 



dramatically in their outer regions, beyond ~ 4 disk scale 
lengths; and (2) at these large distances, the vertical velocity 
dispersions are nearly independent of radius (ctj ~ 20 kms"') 
rather than decreasing exponentially as expected for a con- 
stant mass-to-light ratio, constant scale-height exponential 
disk. These findings are in very good agreement with the the- 
oretical predictions presented in Figures |4] and |5l suggesting 
that the flaring and higher than expected values of cT; in the 
disks of M83 and M94 could be due to bombardment by halo 
substructure. Exploring these constraints further would re- 
quire extending the number of kinematic surveys in the outer 
disks of spiral galaxies as well as performing an extensive se- 
ries of numerical experiments to fully sample the statistical 
variation in halo accretion histories predicted by ACDM. 

In the context of distinctive signatures induced in stellar 
disks by halo substructure, another intriguing result is re- 
ported in Figure|6l This figure illustrates that the face-on sur- 
face density distribution of the simulated disk becomes shal- 
lower at large projected radii (R > 5Rd) as a result of the satel- 
lite impacts, so that there is an excess of light relative to the 
exponential profile of the inner disk. Such behavior is interest- 
ing in the context of the so-called antitruncated disks whose 
surface brightness profiles display similar excesses at large 
dista nces, beyond 4 - 6 disk scale lengths (e .g., Erwin et aU 
lOO llPohlen & Trujillol2006tlPohlen et al.l2007,:,Erwin et all 



20081) . Transport of angular momentum during the subhalo 
bombardment causes the migration of disk material outward 
in radius, leading to the excess surface density at large radii. 
Of course, satellite-disk encounters of the kind considered 
here do not constitute the only mechanism for such angular 
m omentum transfer. 

lYounper et alj (12007 ^ recently investigated the origin of an- 
titruncated disks in face-on spirals using hydrodynamical sim- 
ulations of minor mergers of galaxies. These authors sug- 
gested that the antitrucation is produced by the action of two 
competing effects; (1) merger-driven gas inflows deepening 
the central potential of the primary galaxy and contracting the 
inner surface density profile of the disk; and (2) angular mo- 
mentum and stellar mass transfer in the outer regions, caus- 
ing the disk to expand. The first process acts to maintain an 
inner scale length similar to that of the initial disk of the pri- 
mary, while the second is responsible for generating the ex- 
cess surface density relative to the exponential profile of the 
inner disk. This latter process is similar to what occurs in our 
n umerical experiments. 

lYounger et al.l (2007) further argued that the combination 
of these two mechanisms produces antitruncated disks only in 
the hydrodynamical regime. They argued that in collisionless 
interactions in which the central potential is not as deep, the 
disk expands more uniformly in response to the angular mo- 
mentum transfer. As a result, the entire surface density distri- 
bution becomes flatter and the disk scale length increases at all 
radii. Our experiments are dissipationless, but the deep cen- 
tral potential needed to maintain the inner scale length may 
be provided by the massive, central bulge. Indeed, fitting the 
mass surface density of the final disk to an exponential pro- 
file yields a scale length only 8% larger than that of the ini- 
tial disk. The implication is that the satellite-disk interactions 
that should be common during the hierarchical assembly of 
galaxies may have a non-negligible role in the formation of 
at least some of the antitruncated disks, even when encoun- 
ters are completely dissipationless. Of course, the inclusion 
of gasdynamics, star formation, and stellar components in the 
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infalling systems as well as a detailed comparison of the fi- 
na l disk with ob s erved antitruncated systems, as was done 
in lYounger et alj (120071) . will be important to fully explore 
the origin of antitruncated disks in the context of hierarchical 
CDM. 

Fourier analysis of their surface density illustrates that the 
simulated stellar disks also exhibit lopsidedness spawned 
by the interactions with CDM substructure (Figure |7]l. 
This has important implications as a significant fraction of 
observed spiral galaxies display lopsided asymmetries in 
their gaseous and/or stellar distributions (e.g., Baldwin et d] 
19801: iRic hter & Sancisi 1994; Rix & Zaritskv 1995; 



Zaritskv & Rix 
1998t 



199 7; .Matthews et al. 1998; 
Rudnick & RTxI I l998t ' 

~ 2007t iReichard et ah i200: 



H aynes et ah 
Bournaud et 120051: 
For 



Angiras et a l. 2006, . . . 

example, [Richter & Sancisil (^1994^ analyzed the Hj dis 
1700 disk galaxies and concluded that 
their sample 



tributions m 

the frequency of lopsided systems m tneir sample ex- 
ceeded 50%. Lopsided stellar disks are also ubiquitous 
jRix & Zaritsk vlll995l: IZaritskv & Rix' '1997^, 'Rud nick & Rixl 
[1998; BournaudeTaDllOOSt [Reichai-d et al. 2008^^ Other 
proposed mechanisms for generating lopsidedness have 



included tidal interactions and me rgers (Walker etal. 19961 
Zaritskv & Rixll997tlAngiras et al.i 2006. 2007; Mapelli et^ 



20081). asymmetric accret ion of intergalactic gas into the 



disk dBournaud et al.ll2005h. offsets betw een the disk and the 
dark matter halo (^Levin e & Sparkd '1998). dynamical insta- 
bilities/processes internal to th e disk (Sellwood & Merritt 
1994; Sver & Tremaine 1996; Sellwood & Valluri 1997; 
De Riicke & De battista 2004; Saha et aL _2007: Dury et al.. 
2008h . and halo lopsidedness dJog 1997, il"999ir it is plausible 



that a number of mechanisms operate in concert to foster 
lopsidedness. 

Further comparison of our results with observational work 
on stellar lopsidedness is worthwhile. Zaritsky & Rix (19971) 
used near-infrared photometry to study a sample of 60 late- 
type spiral galaxies in the field. These authors defined a quan- 
titative measure of lopsidedness as the radially-averaged ratio 
of the m = I to OT = Fourier amplitudes between 1 .5 and 
2.5 disk scale lengths Rij, and denoted this quantity by (Ai). 
Averaging over a range of radii reduces the effect of isolated 
asymmetric peaks on the results of the observational analy- 
sis and provides a globa l meas ure of the lopsidedness in each 
galaxv. IZaritskv & Rixl (Il997h found that ~ 30% of all disk 
galaxies in their sam ple exhibit significant stellar lopsided- 
ness with (Al) > 0.2. iRudnick & ^ (Il998h followed up on 
this work by using R-band photometry to investigate asym- 
metries in 54 early-type spirals. These authors reported a me- 
dian value of (Al) in their sample of 0.11 and that 20% of 
their disk galaxies had (Ai) > 0.19. The similar amplitude of 
lopsidedness in disks with very different star formation rates 
indicates that the majority of the observable asymmetries in 
the stellar light of galactic disks reflects asymmetries in the 
stellar mass distribution rather than asymmetric star forma- 
tion, confirmin g a dynamical origin o f stellar lopsidedness. 
More recently, Bournaud et al.l (120051) analyzed a sample of 
149 galaxies from the Ohio State University Bright Galaxy 
Survey observed in the near-infrared and reported a mean (A i ) 
equal to 0.1 1. 

We have followed these authors and computed the quantity 
(Ai) in our simulated stellar disks. All measurements were 
performed after allowing the disks to relax from the encounter 
with each infalling subhalo as described in § 12.41 As shown in 
Figure |6] the slope of the surface density profile and hence the 



disk scale length does not evolve significantly in the relevant 
radial range for measuring (Ai), l.5Rd S '^■5Rd- There- 
fore, we adopted the scale length of the initial disk for these 
calculations and computed the quantity (Ai) in each case by 
averaging A 1 between ~ 4.2 kpc and ^ 7 kpc. We also con- 
firmed that using different bin numbers in calculating (A i ) 
gives similar results. 

For the simulated accretion history of host halo Gi, we de- 
rived (Al) values in the range 0.10< (Ai) <0.16. For ref- 
erence, the initial axisymmetric disk had (Ai) ~ 0.003. We 
stress that much larger values of (Ai) are excited during the 
subhalo impacts where tidal forces are strongest, but these 
are shorter-lived. As expected, the amplitude of (Ai) de- 
creases steadily after the last satellite passage, reaching a 
value of ~ 0.013 in the final disk. We also found that the 
magnitude of the induced lopsided asymmetries depends sen- 
sitively on the structural properties of the infalling subhalos 
as well as those of the galactic disks. For example, satel- 
lites modeled with steep density profiles generated a signif- 
icantly stronger lopsidedness compared to that of the cored 
counterparts. In addition, bulgeless disks develop more pro- 
nounced lopsided asymmetries in their inner regions com- 
pared to disks with massive bulges. These trends can be read- 
ily understood as a consequence of the strength of subhalo 
tidal perturbations onto galactic disks. Extending the mass 
spectrum of infalling subhalos to include the mos t massive 
systems with Mjub > A/disk (e.g., iPurceU et all l20()8l) would 
be required to examine whether satelUte bombardment can 
explain the largest lopsided asymmetries observed in some 
galaxies. Additional mechanisms for lopsidedness, as dis- 
cussed above, may also need to be invoked in this respect. 

The above analysis demonstrates that encounters with 
CDM substructure can excite lopsidedness in stellar disks 
at levels similar to those measured in observational sam- 
jgles of disk galaxies ( Rix & Zaritskv 1995; Zaritskv & Rix! 
1997; Rudnick & Rix 1998). It also indicates that ACDM- 
motivated subhalo accretion histories can maintain these lop- 
sided asymmetries for a significant fraction of the cosmic time 
(> 5 Gyr). This becomes particularly relevant in light of 
the fact that several studies find no observed correlation be- 
tween the p resence of nearby com panions a nd disk lopsid- 
edness(e.g., IZaritskv & Rixl Il997t [Wilcots & PrescottI 12004 
iBournaud et al.ll2005b . 

Our simulations also suggest that significant disk tilting 
may result in response to encounters with CDM subhalos 
(Figure[8]). Disk tilting has potentially important implications. 
The formation of disk galaxies remains poorly understood de- 
spite som^recentadvances^ Weil et al. 1998; Abadi et all 
120031; ISo mmer-Larsen et alj |2003£ .Governato et al.. ,2004t 
Robertson et al.1 12004 iGovernato et al.1 120071) . The relation 
between the angular momenta of galactic disks and the net 
angular momenta of their host halos remains unclear Dark 
matter halos typically have their net angular momenta aligned 
with the shortest of their principle axes (though the degree 
depends on halo mass, e.g.. lBailin & Steinmetzll2005l) . In the 
Local Group, the disks of the MW and M31 seem to have an- 
gular momenta that point along the larger-scale s tructure de- 
lineated by the Local Group dwarf galaxies (Maiewski"1994; 
Hartwick 2000 ; Will man et al. 2004; Zentneretal. 2005a 
ILibeskind et al.l l2005h . Several studies have attempted to 
test for such aUgnments between disk principle axes and 
large-scale structur e (Zaritskv et al .1 1 1 9971: lYang et al.l 120061: 
lAzzaroet all 120061 l2007c iBaiUn et aLn2008l) . but little evi- 
dence for either alignment have been found in statistically- 



CDM Substructure and Galactic Disks II: 



21 



large samples. Our results indicate that in addition to intrinsic 
scatter among the an gular momenta of halos an d their large- 
scale environments ( iBailin & Steinmet3 l2005h such align- 
ment may be diluted because disks, once formed, may tilt in 
response to numerous interactions with infalling satellites. 

As far as the robustness of galactic disks to encounters 
with satellites is concerned. Figures H] |5] and |9] illustrate that 
thicker disks are relatively more resilient to subhalo bom- 
bardment compared to their thin counterparts (see, however, 
ISellwood et alJll998h . Infalling satellites deposit energy into 
galactic disks via impulsively shocking individual orbits of 
disk stars during their passage ("direct heating") as well as 
by exciting global collective modes in the disk. Collective 
modes include both vertical bending waves (e.g., warps) and 
horizontal density waves (e.g., spiral structure and bars) and 
we stress that there is coupling between planar and vertical 
modes in 3D. In this case, the energy is transferred to the disk 
causing its heating by damping of the waves via resonant cou- 
pling (e.g., Weinberg 1991). Alternatively, the waves may be 
damped by dynamical friction exerted by the dark matter halo, 
thus yielding no heating of the disk. 

Calculations of direct heating show that during a single, 
random orientation passage of a satellite with mass Msat mov- 
ing at relative speed v and impact parameter b, the mean 
change in the vertical energy per unit mass of a disk star is 
given by AE, = (hj/3)(GM,Jb^K,fp^L(j3) (lSpitzerill958l) . 
Here, h. is the rms thickness of the disk, is the frequency 
of vertical oscillations, the parameter (3 = In-b/v is of order 
the characteristic passage time of the satellite divided by the 
orbital period of the star, and Lip) is a dimensionless function 
which is unity for /3 ^ and exponentially small for (3':^ I. 
According to this formula, direct heating should vanish as the 
disk becomes razor thin, both because h, and also be- 
cause oo. This fact in conjunction with the results in 
Figures |4]|5] and |9] suggests that the relative fragility of thin- 
ner disks to encounters with satellites lies in both the ability 
of these systems to support global instabilities and the effec- 
tiveness of collective modes to transfer energy in this case. 
The absence of a bar and the much weaker warp and spiral 
structure induced by the subhalo impacts in galaxy model D2 
compared to its thinner counterpart Dl lends support to this 
conjecture. Definitive investigations of these issues would re- 
quire a combination of targeted numerical ex periments as well 
as exten sions of earlier analytic work (e.g.. IWeinberdlT989t 
[Weinbe rg & BHtz 2006). 

A relevant issue concerns the existence of bulgeless, thin 
disk galaxies in cosmological models where accretions of 
massive satellites are as common as predicted in ACDM. Fig- 
ure [To] demonstrates that a bulge component reduces signif- 
icantly the damage done to the disk, so that bulgeless disk 
galaxies experience substantially more thickening by infalling 
satellites compared to their counterparts with bulges. Because 
most MW-sized halos are expected to have accreted numerous 
substructures that are a significant fraction of the disk mass in- 
cluding at least one system as massive as the disk since z = 1 
(Figure [T] Stewart et alj|2008b . the ubiquity of very thin, bul- 
geless disk galaxies containing dominant old stellar popula- 
tions would be difficult to reconcile with ACDM. Interest- 
ingly, using SDSS data K autsch et al. (200 6) recently com- 
piled a uniform catalogue of 3169 edge-on disk galaxies and 
found that ^ 1/3 of the galaxies in their sample were bul- 
geless, "super-thin" disks with extreme axial ratios. More- 
over, all systems in the sample of bulgeless, edge-on spi- 



rals of iDalcanton & BernsteinI ( 120021) have pronounced thick 
disks and there are no signs of companions in the vicinity of 
the prototype thin bulgeless disk galaxy M33. Formulating 
a comprehensive model for the formation and survivability of 
very thin, bulgeless disk galaxies in the context of hierarchical 
CDM remains challenging. 

The findings of the present study as well as those of 
Paper I have interesting implications for the formation 
of thick disks. Thick disks are structurally, chemically, 
and kinematically distinct from thin disks, and there is 
evidence that they ma y ass emble quite early in the history 
of a galaxy dElmegreen & ElmegreenI l2006l) . Thick-disk 
stars in the MW and external galaxies are characterized 
by much larger scale heights, exhibit larger velocity dis- 
persions and slower rotation, and are more metal-poor 
and significantly en hanced in g-elements compared to 
thin-d is k stars (e .g., Reid & Maiewski 1993; Gilmore et alJ 
1^951 iWvse & Gilmore 1995; Prochaska et al. 200a 
Chiba& Beers 2000; Bensbv et al. 2005; Seth et al. 2005; 
Yoachim & Dalcanton 2006; Al lende Prieto et al., ,20061; 
Juric et al. 2008; Ivezic et al. 2008). While our dissipation- 
less simulations can neither verify nor disprove any of the 
trends regarding metallicities, the present study does show 
that encounters with infalling subhalos increase considerably 
the scale heights and velocity ellipsoids of thin, galactic disks 
(Figures]?] and |5]l. Furthermore, the vertical structure of the 
final disk is well-described by a standard "thin-thick" disk 
decomposition (see Figure ]4] in Paper I) and analysis of the 
mean azimuthal velocity of disk stars at the solar radius in our 
simulations also reveals a vertical gradient in rotational veloc- 
ity of ^ -20 kms"^ kpc~' between 1 and 3 kpc from the disk 
plane, which is consistent with what is inferred for the thick 
disk of the MW (Allende Prieto et al. 2006). These results 
suggest that at least part of a galaxy's thick-disk component 
may plausibly originate from the vertical dynamical heating 
of preexisting thin disks by CDM substructure. While this 
conclusion is supported by a number of o bservational studies 
in both the MW ( e.g., ]Robin et al.1 [1996; Chenetal. 200|t 
iBensbv etal] 1200 5') and external galaxies (iSeth et al.ll2005l) . 
more detailed theoretical modeling of the properties of thick 
disks would at least require the inclusion of gasdynamics, 
star formation, and metal enrichment. 

Of course, vertical dynamical heating of an existing thin 
disk does not constitute the only viable model for the ori- 
gin of thick disks. Other proposed mechanisms include 
satellite accretion events that directly deposit thick-disk stars 
at large scale heights (e.g. .Statler 1988; Abadi et al. .2003]; 
I Yoachim & Dalcantonl2005l 120081) as well as thick disks stars 
forming in situ at early times dire ctly from gas at la r ge dis - 
tances above the midplane (e.g.. Brook et al] 12004 l2005h . 
While there is no definitive observational or theoretical evi- 
dence to rule out any of the models for the origin of thick disks 
conclusively, the existence of very slowly rotating or counter- 
rotating thick disks in a significant fraction of disk galaxies 
would be pa rticularly problematic for the "vertical heating" 
mechanism (Yoachi m & DalcantonI |20()5]) . Most likely, all 
mechanisms do operate simultaneously at a different degree 
in forming the thick disk components of galaxies. 

In this paper, we have focused on the gravitational inter- 
action between galactic disks and CDM substructure in the 
collisionless regime. Given the complex interplay of effects 
relevant to the formation and evolution of galactic disks, the 
inclusion of gas dynamics, star formation, and chemical evo- 
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Fig. 13. — Effect of the adiabatic growth of a thin-disk component on the vertical structure of a thick disk. For the latter, we adopt disk model D2. Solid 
lines show results for the initial uncontracted thick disk with mass M^i^k while the remaining lines correspond to the final state of the thick disk after the growth 
of thin disks with various masses, Mg,ow All growing thin disks have a sech^ scale height of zthin = 0.2 kpc and the same radial scale length, Rj, as that of 
the initial thick disk. Left: Disk thickness profiles. Thicknesses and radii are normalized to the scale height, zd, and radial scale length, Rj, of the initial thick 
disk. Right: Vertical velocity dispersion profiles normalized to the total circular velocity of disk model D2 at the solar radius, Vc(Rq) = 235.6 kms"'. The slow 
accumulation of thin disk material modifies the structure of the initial thick disk, causing its vertical and radial contraction as well as an increase of its vertical 
velocity dispersion. 



lution would be required at the least to refine the conclusions 
presented here. Specifically, the modeling of hydrodynamics 
will be crucial in determining the extent to which the presence 
of gas can influence the dynamical response of galactic disks 
to satellite accretion events and affect the properties of the fi- 
nal disk. For example, both subhalos and the disk may contain 
gas, particularly at high redshift. The satellite accretion events 
would then trigger bursts of star formation that may replenish 
a thin disk component. 

Moreover, a dissipative component may alter the dynam- 
ical effects of substructure on stellar disks in two important 
ways. First, the gas itself can absorb part of the orbital en- 
ergy deposited into the galactic disk by the infalling systems, 
acting as an energy sink. This process would lead to the heat- 
ing of the gaseous component. However, the efficiency of this 
mechanism in reducing the dynamical damage done to a stel- 
lar disk will depend on the gas content of galactic disks at 
early times when most of these accretion events occur. Inter- 
estingly, analytical models for the evolution of the MW disk 
in a cosmological context do estimate that the gaseous disk at 
I s hould amount to ^ 50% of the mass of the stellar disk 
jNaab & Ostriker 2006). Given that substantial gas fractions 
are expected at high redshifts, the role of gas in stabilizing 
the galactic disks against the violent gravita tional encounters 
with satellites may be crucial (IStewart et al.l l2009). 

Second, owing to its dissipative nature, the gas can radiate 
its energy away. As a result, the gas that has been heated by an 
encounter with a subhalo can subsequently cool and reform a 
thin disk. As any gaseous component slowly accumulates in 
the center of the mass distribution, it will also induce con- 
comitant contraction of the thickened stellar disk due to its 
gravity. Larger s cale smooth gas ac cretion in galaxy forma- 
tion models (e.g.. iMurah et al.ll2002l) will also be relevant in 



this context. A full exploration of these contingencies is chal- 
lenging and we defer such studies to future work. In what 
follows, we present a simple, preliminary experiment to serve 
as a crude estimate of such effects. 

In particular, we investigated the response of an initial thick 
disk to the adiabatic growth of a massive, thin-disk compo- 
nent within it. For the former, we adopted galaxy model D2 
with a vertical scale height of Zd = I kpc. We considered 
three growing, exponential thin disks with masses Mc,iow = 
[0.1,0.5,1 ]Mdisk, where Mdisk is the mass of the disk in model 
D2, and a sech" scale height of Zthin = 0.2 kpc. The latter 
value is consistent with the scale heights of known, young, 
star-forming disks observed in both e xternal galaxies (e.g., 
'Wainscoat et al."'1989';'Matthews"2000^ and in the MW (e.g., 
Bahcall & Soneira 1980; Reid & Majewski 1993). AH grow- 
ing disks were treated as rigid potentials and had the same 
radial scale lengths as disk D2. Furthermore, the scale length 
of each thin disk was kept constant while its mass was slowly 
increased from zero to its final value linearly over a timescale 
of 1 Gyr. Such timescale s are in general accordan ce with disk 
formation models (e.g.. F all & Efstathioulll980l) and ensure 
that the process of disk growth is approximately adiabatic. 

Figure [13] presents the results of these experiments. In all 
cases, the initial thick disk contracts vertically as well as ra- 
dially in response to the growth of the thin-disk component. 
Moreover, its vertical velocity dispersion increases, reflecting 
the deepening of the potential well due to the slow accumu- 
lation of thin-disk material. As expected, these changes in 
the structure of the thick disk depend sensitively on the total 
mass of the growing disk. In the most dramatic case with 
Mgiow = A^disk, the decrease in the thickness of disk model 
D2 is ^ 30% in the solar neighborhood. Furthermore, these 
changes do not occur uniformly as a function of radius. The 
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disk potentials are centrally concentrated, so the evolution of 
the inner disk is much more pronounced compared to that of 
the outer disk. Overall, the results presented in Figure[T3]sug- 
gest that the mass of the growing disk would need to be many 
times that of the thickened stellar disk to reduce its thickness 
appreciably, and even then, there would be stability issues. 

It is important to emphasize that the aforementioned ex- 
periments were not designed to elucidate the importance of 
adiabatic thick-disk contraction in the specific case of the 
MW whose midplane density ratio of thick-to-thin disk is 
only 12% (e.g., Juric et al. 2008). In a study of thick disks 
in the Hubble Space Te lescope Ultra Deep Field (UDF), 
lElmegreen & ElmegreenI (i2006) did examine whether such a 
process could have determined the present-day scale height 
of the thick disk of the MW. Their calculations showed that 
if the present thick-disk component of the MW began as an 
equilibrium pure-thick disk at a young age, and if subsequent 
accretion of the entire thin disk was adiabatic, then the initial 
thick-disk scale height had to be ^ 3 kpc. This is consid- 
erably larger than that observed for young thick disks in the 
UDF, where the average scale height is 1 .0 ± 0.4 kpc. 

6. COMPARISON TO PREVIOUS WORK 

The response of disks to encounter with infalling satellites 
has received a great deal of attention owing to the numerous 
implications that it entails for theories of galaxy formation 
and evolution. Here we discuss the main differences between 
o ur results and thos e re ported in a subset of p revious studies. 

iFont et alJ (120011) and lGauthier et al.l C2006) carried out nu- 
merical simulations of the gravitational interaction between 
galactic disks and a large ensemble of dark matter subhalos. 
Both of these investigations considered the z = Q satellite pop- 
ulations present in a MW-sized CDM halo, and both stud- 
ies reached the conclusion that halo substructure has an in- 
significant ef fect on the global structure of a galactic disk. 
In particular, iFont et af.l (12001 ) showed similar tidal heating 
rates in the ir A^-body stell a r disk with and without substruc- 
ture, while iGauthier et al.l (l2006h reported appreciable heat- 
ing in the inner disk regions due to the formation of a bar and 
only mild heating of the disk at intermediate and large radii. 
These results are explained b y the fa ct that the strategies of 
[Font et al.1 (|2001) and Gauthieretal] (|2006) excepted those 
systems that are most capable of strongly perturbing the disk. 
Massive subhalos on highly eccentric orbits at early epochs 
suffer substantial mass loss or become preferentially disrupted 
during their orbital evolution in the host potential p rior to 
z = (Zentner & Bullock 2003; Gao et al. 2004; Zentne r et all 
l2005al ; iBensonl i2005), and so they are more likely to be ab- 
sent from t he pre sent-day subhalo populations. Indeed, in the 
iFont et akl (1200 Ih experiments only subhalos with masses be- 
low 1O^M0 had pericenters within the solar radius, Rq. A 
primary improvement in our modeling is that we have fol- 
lowed the formation history of the host halo since z ^ I, and 
consequently we have accounted for a larger number of im- 
portant satellite-disk interactions than that based on the z = 
substructure. As a result, we report significantly more dam- 
age to the structure of the galac tic disk than that de monstrated 
by either F ont et al.l (2001) or Gaut hier et aP (2006). 

Past numerical investigations into the resilience of galactic 
disks to infalling satellites have also suffered from nu- 
merical limitations and/or from assumptions which curtail 
their ability to accurately capture the degree of global 
dynamical evolution that accreting subhalos can induce in 
cold, stellar disks in a cosmological context. For exam- 



ple, the modeling of various components in the primary 
disk galaxy and/or the s atellites as rigid potentials (e.g., 
'Ouinn & Goodman" 1986"; 'Ouinn etal."1993'; 'Sell wood et al] 
1998; Ardietal. 2003; Hayashi & Chiba 2006), the ini- 
tialization of a disk much thicker compared to typical thin 
disks such as the old, thin s tellar disk of t he MW (e.g., 
Ouinn & Goodman 1986; O uinn et alj 119931; IW alker et al] 
1996; Huang & Carlberg 1997t I Velazquez & White , 19991 
Fontet al. 2001; Villalob os & Helmill2008l) .' and the infall of 
satellites with orbital parameters and/or structural properties 
inconsistent with ACDM expectations (e.g., Ouinn et aU 
Walker et alj__p96t OHuang & Carlberg 1991 
'Velazquez & White 1999|). For ex ample. Ouinn & GoodmanI 
(1986), Ouinn etal. (19911), and IWalker etalJ ( Il996h per^ 
formed the first numerical explorations of the interaction 
of single satellites with larger disk galaxies. These studies 
unanimously found disks to be quite fragile, reporting that 
an encounter with a satellite of only 10% of the disk mass 
could increa se the disk thickness by a factor of -^^ 2 at the 
solar radius dOuinn & Goodmanlll986l;IOuinn et"anil993h . In 
contrast, we find galactic disks to be generally more robust 
to accretion events than these earlier investigations have 
indicated. 

These differences may be due to a variety of factors. When 
a satellite is modeled as a distribution of interacting particles 
as opposed to a rigid potential, the efficiency with which it 
can heat a galactic disk is suppressed for two reasons. First, 
a self-gravitating satellite suffers mass loss due to tidal strip- 
ping and shocking. Second, a responsive satellite can absorb 
part of its orbital energy, decreasing the amount of energy de- 
posited into random motions of disk stars. Moreover, live ha- 
los are needed to treat the effect of an accreting satellite on 
disk structure properly. Representing DM halos as rigid po- 
tentials leads to an over-estimate of disk thickening by a fac- 
tor of ~ 1.5-2, whereas a self-gravitating halo will respond 
to both the disk and satel l ite and aid in stab ilizing the disk 
dNelson & Tremaindl 1995 1 ; I Velazquez & Whit&.1999) . More- 
over, the focus on prograde circular or nearly circular orbits 
in some of the aforementioned studies is also likely to have 
overestimated the typical amount of disk damage. The orbit 
with the most damaging effect for a thin disk is a coplanar, 
prograde circular orbit, since this causes the satellite force to 
be in near resonance with the disk stars. More eccentric orbits 
are likely to cause less damage to a disk. This is immediately 
apparent from the fact that, in the impulsive limit, the energy 
transfer scales as oc v"^. To a first order, the high-speed peri- 
centric passages of CDM subhalos on highly eccentric orbits 
are thus expected to cause less thickening than circular ones. 
However, it is unclear to what extent highly eccentric orbits 
can excite global modes in the disk. This requires a detailed 
series of numerical experiments to evaluate. 

Velazquez & White ( 1999) conducted a large number of A^- 
body experiments in which they investigated the dependence 
of disk heating on, among other variables, the inclination and 
eccentricity of the initial satellite orbit. These authors showed 
that prograde encounters are more efficient at thickening the 
disk than retrograde ones, arguing for a resonant coupling 
between the satellite and disk ang ular momenta th at is sup- 
pressed in retrograde interactions. Velazquez & White (19991) 
found that a satellite with an initial mass of Msub = 0.2Mdisk 
on a prograde orbit causes an additional increase in the scale 
height of the initial disk at the solar radius of ^ 35% com- 
pared to its retrograde counterpart. 
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While our simulations indicate qualitatively similar dif- 
ferences between prograde and retrograde interactions (Fig- 
ure [T2] |. they show a less sig n ifican t differential effect com- 
pared to IVelazquez & Whita (Il999h . This discrepancy has 
two plausible sources. First, we study cosmologically- 
motivated, comparably high-speed passages of substructures. 
In such rapid ("impulsive") interactions, the resonant cou- 
pling between the satellite and disk angular momenta is ex- 
pected to be weaker, explaining why the relative effect of pro- 
grade versus retrograde interactions is smaller in our exper- 
iments. Second, we study single subhalo passages, whereas 
IVelazquez & White! (1999) followed the infalling satellites 
over many orbits until they merge with the galactic disk or 
become disrupted in the p rocess. This is relevant because 
IVelazquez & Whit3 (119991) found that during the early stages 
of the interaction when the satellite orbit is still eccentric the 
difference between the amount of disk thickening caused by 
prograde and retrograde encounters was fairly small. Satel- 
lites on prograde orbits constituted much more efficient per- 
turbers compared to their retrograde counterparts only during 
the late stages of the interaction when, owing to dynamical 
friction, the satellite speed has been reduced, and its orbit has 
both circularized and become coplanar with the disk (H. Ve- 
la zquez 2008, priv ate comm un ication). 

iToth & Ostrikeri ([1992h . iBenson et al.] (120041) . and 
[Hopkins et al.l (12008 ) quantified the fragility of galactic 
disks to infalling satellites by using semi-analytic models of 
different degrees of sophistication. Analytical approaches 
have the advantage of not being limited by numerical resolu- 
tion, allowing the calculation of a statistically-large number 
of model realizations, but also have the drawback that they 
cannot account fully for the non-linear interaction between 
satellites and disks. Yet, as both our results in Figure [TT] 
suggest and other more targeted studies have indicated (e.g., 
[Weinberg 1991; Sellwood et al. 1998), the excitation of 
global collective modes in the disk is an essential mechanism 
f or energy deposition by a n accreting satellite. 

iToth & Ostrikerl (Il992h concluded that sinking satellites 
with only a few percent of the disk mass could lead to sub- 
sta ntial disk thicken ing. Specific comparisons with the work 
of iToth & OstrikeJ ([1992) is very difficult. In a sense, the 
fact that their disk suffers substantial damage strengthens our 
findings. However, their results are averaged all possible or- 
bital inclinations making specific model comparisons cumber- 
some. Another difficulty lies in their assumption that the or- 
bital energy of the satellite is deposited locally at the point of 
impact. Our results do indicate that the energy imparted by 
typical cosmological substruc tures will be deposited globally 
across the entire disk. Lastlv. .Toth & Ostrikerl (1 19921) assumed 
that the total thickening and heating scale with satellite mass 
and that they are the same regardless of the initial thickness 
and velocity ellipsoid. Figures |4] and |9] show that this con- 
jecture was also incorrect. The implications of most of these 
assumptions ar e unclear. 

In contrast, iBenson et alj (120041) suggested that the ob- 
served thickness of stellar disks is entirely compatible with the 
abundance of substructure in CDM hal os. In addition to no t 
accounting for global collective modes, IBenson et al.l (120041) 
also adopted satellite orbits that spanned only a limited range 
of orbital energies and angular momenta in spherical halo po- 
tentials rather than the richer variety of impacts e xperienced 
over the cours e of halo formation (e.g.. lZentner et al. 2005a b; 
lBensonir2005h . Consequently, this study suffers from a flaw 
similar to that of .Font et al.. (.2001.) and .Gauthier et al.. (.2006h . 



explaining possibly why we report more significant damage 
to the disk structure by subhalo bombardment. 

More direct comparison can be performed with the recent 
work of Hopkins et al. (2008). These authors have argued 
that deposition of orbital energy into the disk in the context 
of realistically radial subhalo orbits and ACDM-motivated 
accretion histories would yield a much less dramatic im- 
pact on the disk structure than previously th ought, with or 
withou t the presence of gas. In particular, iHopkins et alj 
(l2008h derived that the disk heating efficiency is nonlinear 
in mass ratio, cx (Mjub/A^disk)^, instead of the linear scal- 
ing of Toth & Ostriker (1992), implying that the fractional 
change in disk scale height should be very small even for 
the very massive subhalos w e considered in the pre sent study 
(0.2Mdisk < Msub < Mdisk). 'Hopkinsetal] dlOOSl) defined a 
disk thickening parameter AH/R, where H is the median disk 
scale height and R is the radius where the scale height is mea- 
sured, typically within a factor of two of the disk half-mass 
radius, R/j. 

For the sim ulated accretion history of host halo Gi, the 
IHopkins et all (|2008) formula predicts AH/R ~ 0.01, if we 
simply add the masses of all cosmological subhalos S1-S6. 
Our numerical experiments indicate significantly more thick- 
ening. Indeed, at R = Rh= 1 JRd, we measure AH/R ~ 0.02, 
and because of the pronounced flaring in response to the 
accretion events, we report even larger AH /R ~ 0.04, at 
R = 2Ri,. A variety of reasons may be responsible for these 
discrepancies. As before, global collective modes which 
may dominate the disk response to accretion events were 
not included in these simple analytic scalings. Furthermore, 
IHopkins eFa l. ( 2008 ) calibrated their r esults to numerical 
simul ations d Velazquez & White! 119991: IVillalobos & Helmi! 
12008!) with initial disks that were significantly thicker com- 
pared to the thin, galactic disk we employed here, and were 
therefore intrinsically more robust to encounters with satel- 
lites (Figure fTTTi. 

Lastly, special emphasis s hould be placed on the recent 
study bv lPurcell et al.! (120081) . These authors performed col- 
lisionless A^-body simulations to study the response of an 
initially-thin (zti ~ 400 pc), MW type disk galaxy to 1 : 10 
satellite impacts. Such accretion events represent the primary 
concern for disk survival in the ACDM cosmological model 
and should have been conimonp la ce in the history of G alaxy- 
sized halos dStewart et alj|2008!) . iPurcell et all (12008") quan- 
titatively demonstrated for the first time the destruction of 
a thin, stellar disk by these cosmologically motivated com- 
mon events. They find that regardless of orbital configura- 
tion, the impacts transform the disks into structures that are 
roughly three times as thick and more than twice as kinemat- 
ically hot as the observed dominant old thin disk component 
of the MW. On the other hand, our work shows that a thin 
disk component may survive, even strongly perturbed, the en- 
counters with halo substructure (Paper I). This is because we 
have focused on infalling systems with masses in the range 
0.2Mdi,sk ^ A^sub ^ A^disk, ignoring the most massive accretion 
events that could prove ruinous to thin-disk survival. Our find- 
ings are thus relevant to systems that have already experienced 
the most destructive events simulated by Pur cell et alj (120081) 
since z = 1 and have re-grown their thin disks since. In this 
sense, our simulation set a nd results are complementary to 
those of IPurcell et all dlOOSh . 
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7. SUMMARY AND FUTURE DIRECTIONS 

Using a suite of high-resolution, fully self-consistent dis- 
sipationless A^-body simulations we have examined the dy- 
namical effects of halo substructure on thin galactic disks in 
the context of the ACDM paradigm of structure formation. 
Our simulation campaign utilizes cosmological simulations 
of the formation of Galaxy-sized CDM halos to derive the 
properties of infalling subhalo populations and controlled nu- 
merical experiments of consecutive satellite encounters with 
an initially-thin, fully-formed disk galaxy. As a corollary, we 
have quantified the importance of various physical effects that 
could influence the response of a galactic disk to substructure 
accretion events. The properties we address are the initial disk 
thickness, the presence of a bulge component in the primary 
galaxy, the internal density distribution of the infalling sys- 
tems, and the relative orientation of disk and satellite angular 
momenta. 

Our work expands upon past numerical investigations into 
the dynamical response of galactic disks to merging satellites 
in several ways. One improvement concerns the more realistic 
treatment of the infalling subhalo populations. Previous stud- 
ies of the interaction between disks and ensembles of subhalos 
considered onl y the z = Q surviving substructure present in a 
CDM halo (Fo nt et alJl200ltlGaufliier etani2006l) . This leads 
to estimates of the damage done to the disk that are biased low, 
because massive subhalos with small orbital pericenters that 
are most capable of strongly perturbing the disk are preferen- 
tial removed from the satellite populations over time. We have 
accounted for satellite-disk interactions with merging subha- 
los that typically do not survive to the present day but never- 
theless cause significant damage to the disk. This is the major 
conceptual advantage of our work and this difference drives 
our initial disks to be more dramatically affected by halo sub- 
structure than those in earlier studies. Second, the primary 
disk galaxy models we utilize are fully self-consistent parti- 
cle realizations derived from explicit DFs and are designed 
to satisfy a broad range of observational constraints available 
for actual disk galaxies. Finally, we represent satellite sys- 
tems by equilibrium numerical realizations, whose properties 
(mass functions, orbital parameters, internal structures, and 
accretion times) are extracted directly from cosmological sim- 
ulations of the formation of Galaxy-sized CDM halos. 

Our main conclusions can be summarized as follows. 

1 . Close encounters between massive satellites and galac- 
tic disks since z = 1 are common occurrences in the 
ACDM cosmological model. Statistics of four Galaxy- 
sized CDM halos indicate that, on average, ^ 5 sub- 
halos more massive than 0.2Mdi,sk, where Mjisk is the 
present-day mass of the stellar disk of the MW, pass 
within ^ 20 kpc from the centers of their host halos 
in the past ~ 8 Gyr In contrast, very few satellites in 
present-day substructure populations are likely to have 
a significant impact on the structure of a galactic disk. 
This is because massive subhalos on potentially dam- 
aging orbits suffer severe mass loss or become tidally 
disrupted prior to z = as a result of penetrating deeply 
into the central regions of their hosts (§ |2.2| l. 

2. A conservative subset of one host halo accretion his- 
tory was used to seed controlled subhalo-disk encounter 
simulations. The specific merger history involved the 
accretion of six satellites with masses, orbital pericen- 
ters, and tidal radii of 7.4 x lO'' <Msub/MQ < 2 x lO'", 



'"peii S 20 kpc, and rtid > 20 kpc, respectively, since z = 
1 . These events severely perturb an initially-thin (zd = 
0.4 kpc), MW type disk galaxy (Mdisk ~ 3.5 x IO'^'Mq) 
and imprint the following distinctive dynamical signa- 
tures on its structure and kinematics: 

• The development of non-axisymmetric structures 
including a warp, a moderately strong bar, and 
extend ed ri ng-like features in the outskirts of the 
disk (§01. 

• Considerable thickening and heating at all radii, 
with a factor of 2 increase in disk thickness and 
velocity ellipsoid at the solar radius (§ 13.21 and 
§[I3. 

• Prominent flaring associated with an increase of 
disk thickness greater than a factor of 4 in the disk 
outskirts (§[T2]i. 

• Surface density excesses at large radii, beyond ^ 
5 disk scale lengths, as observed in antitruncated 
disks (§[3311. 

• Long-lived, lopsidedness at levels similar to those 
measured in observational samples of disk galax- 
ies (§[33. 

• Substantial tilting of the disk from its initial orien- 
tation in the host halo, resulting in a misalignment 
between hal o an d disk principal axes and angular 
momenta (§ 13.61 ). 

3. Detailed predictions for the dynamical response of 
galactic disks to subhalo bombardment are subject to 
a variety of assumptions including the initial structures 
of the disk and infalling systems, the prominence of the 
bulge component in the primary disk, and the relative 
orientation of disk and satellite angular momenta (§[4|. 

We close with a few words of caution and a discussion of 
fruitful directions for future work that may lead to more con- 
clusive statements about the detailed structure of disk galax- 
ies. We reiterate that we have only addressed the gravitational 
interaction of galactic disks and halo substructure in the col- 
lisionless regime. A full consideration to the rich structure of 
perturbed galactic disks is challenging and would require de- 
tailed knowledge of galaxy star formation histories, gas cool- 
ing and feedback, among other things. However, such studies 
that treat both the gaseous components of the disk and sub- 
halos, and accreted stars will be fundamental in refining our 
understanding of disk galaxy evolution and we plan to extend 
the present investigation in this direction. 

Specifically, spiral galaxies contain atomic and molecular 
gas in their disks which can absorb and subsequently radi- 
ate away some of the orbital energy deposited by the sink- 
ing satellites. Thus, including a dissipative component in 
the primary galaxy will allow for an estimate of the extent 
to which the presence of gas can reduce the damage done to 
disks and stabilize them against the violent gravitational en- 
counters with satellites. Adding star formation as a further in- 
gredient will provide the opportunity to analyze any reformed 
thin disk and establish the degree to which contraction sub- 
sequent to gas cooling can decrease the thickness of a heated 
disk (Figure [T3]l. The effects of gas dynamics and star forma- 
tion will offer the possibility to determine the magnitude of 
starbursts induced in the disk as a result of subhalo bombard- 
ment, while gaining a deeper understanding of the build-up of 
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the inner stellar halos of galaxies dBullock & Johnstonl2005h . 
Lastly, the inclusion of baryonic components in the satellites 
will enable studies whereby, original disk stars, newly formed 
stars, and accreted stars can be located and studied in their 
final configurations. 

Such predictions will be vital as instruments and surveys 
like SDSS III, RAVE, GAIA, SIM, Pan-STARRS, LSST, and 
TMT are poised to provide spatial and kinematic maps of the 
MW and other local volume galaxies to unprecedented detail 
and depth. Our simulations suggest that these experiments 
should uncover detailed disk structure that is substantially per- 
turbed via interactions with infalling satellites. Our ability to 
interpret these data sets will rely on a comprehensive set of 
theoretical predictions for how galactic disks respond to ac- 
cretion events and how this process is convolved with direct 
stellar and gaseous accretions (either via the same encounters 
or other delivery process). In this sense, detailed probes of the 
local volume offer a valuable and unique avenue for constrain- 
ing the process of disk galaxy formation and galaxy formation 
in general. 
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